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ADAPTIVE  METHODS  FOR  COMPRESSIBLE  FLUID  DYNAMICS 
AFOSR-86-0148  FINAL  REPORT 


Marsha  J.  Berger 

Two  major  research  efforts  have  been  supported  by  this  grant.  The  first  is  the  development  of  an 
adaptive  algorithm  for  hyperbolic  conservation  laws  with  simple  physical  geometry.  This  work  is  based  on 
a  combination  of  two  approaches  -  an  adaptive  mesh  refinement  technique  that  concentrates  computational 
effort  where  it  is  most  needed,  and  a  high  order  Godunov  method  developed  for  high  Mach  number 
compressible  flow.  This  cppiuach  has  aided  in  the  resolution  of  the  weak  von  Neumann  paradox  in  shock 
reflection.  It  was  used  to  perform  the  first  calculation  of  Kelvin  Helmholtz  instability  along  the  slip  line  in 
ramp  reflection  off  an  oblique  wedge.  When  combined  with  other  algorithms,  for  example,  multifluid  track¬ 
ing,  it  was  used  to  compute  the  interaction  of  a  supernova  remnant  with  an  interstellar  cloud,  and  to 
categorize  refraction  patterns  when  a  shock  hits  an  oblique  material  interface.  When  combined  with  an 
elliptic  grid  generator,  it  was  used  to  study  the  diffraction  of  a  shock  over  an  obstacle.  In  each  of  these 
cases,  this  approach  to  time-dependent  fluid  flow  yielded  factors  of  10  to  100  improvement  in  efficiency 
over  equivalent  fine  grid  calculations  with  uniform  resolution.  This  work  was  done  in  collaboration  with 
Prof.  Phil  Colella,  at  Berkeley,  and  Dr.  John  Bell  at  Lawrence  Livermore  Laboratory. 

During  the  period  of  this  grant,  the  two  dimensional  algorithm  was  developed  and  a  paper  was  pub¬ 
lished.  The  code  is  currently  being  prepared  for  release  to  other  users.  The  algorithm  was  recently  extended 
to  three  dimensions.  Due  to  the  simple  data  structures  and  the  use  of  nested  grids  in  the  same  topology  as 
the  coarse  grid,  this  involved  little  additional  algorithm  development.  The  most  major  change  involved  the 
development  of  a  new  fine  grid  generator.  Memory  usage  for  three  dimensional  calculations  is  at  a  prem¬ 
ium.  The  old  (two-dimensional)  grid  generator  often  resulted  in  overlapping  subgrids  or  not  very  efficient 
grids  (encompassing  too  many  cells  that  were  unnecessarily  refined).  After  extensive  experimentation,  we 
arc  now  using  an  algorithm  based  on  edge  detection,  borrowed  from  the  computer  vision  and  pattern  recog¬ 
nition  literature,  to  do  a  smart  decomposition  of  the  underresolved  regions  into  efficient  sub-rectangles. 
One  paper  on  this  work  will  soon  appear,  and  another  is  nearing  completion. 

The  second  major  effort  supported  by  this  grant  is  the  development  of  a  Cartesian  grid  method  to 
solve  problems  in  complex  geometry.  This  will  use  all  of  the  machinery  developed  above,  i.c.  the  high 
resolution  Godunov  methods  and  adaptive  mesh  refinement,  along  with  special  difference  schemes  that 
need  to  be  developed  for  the  irregular  cells  along  the  edge  of  the  domain  where  a  grid  intersects  a  solid 
body.  This  work  is  in  collaboration  with  Prof.  Randall  LeVcque. 

Our  first  two  approaches  to  the  "small  cell"  problem  at  the  irregular  cells  used  the  large  time  step 
wave  propagation  method  of  LeVeque,  and  the  flux  redistribution  algorithm  of  Chcm  and  Colella.  Prob¬ 
lems  with  both  of  these  algorithms  (reduced  stability  in  the  first  case,  and  poor  accuracy  in  the  second)  led 
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us  to  develop  a  third  third  approach  to  the  "small  cell"  problem.  We  have  developed  a  rotated  difference 
scheme  for  use  in  the  normal  and  tangential  directions  at  each  point  in  the  boundary,  and  are  currently 
extending  it  to  second  order  accuracy.  In  addition  to  algorithm  development,  this  projects  leads  naturally 
to  research  into  the  accuracy  of  difference  schemes  on  irregular  grids.  This  work  has  relevance  to  other 
types  of  approaches  to  this  problem,  for  example,  unstructured  mesh  algorithms  for  complex  geometries. 
Currently  we  have  a  two  dimensional  implementation  of  a  Cartesian  mesh  algorithm  for  time-dependent 
flow,  including  an  initial  implementation  of  the  adaptive  mesh  refinement  algorithm  described  above. 
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ABSTRACT 

We  describe  a  Cartesian  mesh  algorithm  with  adaptive  mesh  refinement  for  com¬ 
puting  fluid  flows  in  complicated  geometries.  Stable  boundary  conditions  are  needed  at 
the  irregular  cells  where  the  Cartesian  mesh  intersects  the  body.  We  develop  a  difference 
scheme  that  is  stable  even  when  these  irregular  cells  are  orders  of  magnitude  smaller  than 
the  regular  cells.  We  illustrate  its  performance  with  some  computational  examples  solv¬ 
ing  the  two  dimensional  Euler  equations  for  inviscid  flow. 


Introduction 

We  describe  a  Cartesian  mesh  method  to  solve  fluid  flow  problems  in  complicated  geometries.  In 
this  approach,  we  keep  a  uniform  rectangular  (Cartesian)  grid  and  allow  the  solid  boundary  to  intersect  the 
grid  cells  in  an  essentially  arbitrary  way.  Cartesian  meshes  are  an  appealing  way  to  simplify  the  grid  gen¬ 
eration  problem  for  complex  domains.  Multiply  connected  domains  and  irregular  geometries  are  only 
slightly  more  complicated  than  a  simple  domain. 

Cartesian  meshes  have  by  and  large  been  overlooked  in  favor  of  body-fitted  meshes  or  the  more 
recently  popular  unstructured  meshes,  but  they  deserve  much  more  attention.  Cartesian  meshes  have  the 
advantage  of  allowing  the  use  of  high  resolution  methods  for  shock  capturing  that  are  difficult  to  develop 
on  unstructured  grids.  They  also  allow  for  efficient  implementation  on  vector  computers  without  using 
gather-scatter  operations  (except  at  boundary  cells).  They  incur  little  computational  or  memory  overhead 
since  there  are  no  metric  terms  and  they  use  far  fewer  pointer  arrays  than  their  unstructured  counterparts 
Among  the  few  references  on  Cartesian  mesh  methods  are  [Clarke,  Salas  and  Hassan;  Choi  and  Gross- 
man], 

The  major  technical  issue  in  Cartesian  mesh  methods  is  the  small  cell  problem.  Arbitrarily  small 
cells  arise  at  the  edge  of  the  domain  where  the  grid  intersects  a  body.  Stable,  accurate  and  conservative 


difference  schemes  are  needed  for  these  cells.  Moreover,  the  time  step  for  a  time-accurate  computation 
should  still  be  based  on  the  volume  of  the  regular  cells  away  from  the  body,  and  not  restricted  by  small 
boundary  cells.  Previous  efforts  to  use  Cartesian  meshes  have  merged  these  small  cells  together  until  a  cell 
with  sufficient  volume  for  stability  is  obtained.  Clearly  this  loses  resolution.  In  addition,  if  cell  size  is  not 
taken  into  account  in  implementing  finite  difference  schemes  on  irregular  grids,  a  second  order  scheme  can 
lose  one  or  two  orders  of  accuracy. 

Our  treatment  of  boundaries  can  be  combined  naturally  with  our  adaptive  refinement  strategy  using 
locally  uniform  meshes.  We  retain  the  advantages  (efficiency  and  accuracy)  of  uniform  grids  and  are  able 
to  resolve  fine  scale  flow  features  induced  by  complex  geometries.  We  are  using  the  adaptive  mesh 
refinement  algorithm  (AMR)  described  in  [Berger  and  Colella]  to  achieve  accuracy  comparable  to  the 
body-fitted  meshes,  where  grid  points  can  be  bunched  in  an  a  priori  manner  to  improve  the  accuracy  of  the 
solution. 

A  more  complete  description  of  our  approach  to  developing  stable  boundary  conditions  for  irregular 
cells  is  described  in  [Berger  and  LeVeque],  Here,  we  only  sketch  the  main  ideas,  and  present  computa¬ 
tional  examples  for  two  dimensional  time  dependent  flow  in  several  different  geometries. 

One  Dimensional  Model  Problem 

We  motivate  the  basic  approach  in  one  dimension,  where  we  solve  the  equation  u,+f  (u)x  =  0  on  a 
uniform  grid  except  for  one  small  cell  in  the  middle  (see  Figure  1).  Let  h  be  the  cell  size  of  the  uniform 
grid,  and  ah  the  small  cell  size,  0  <  a  <  1.  We  use  an  explicit  finite  volume  scheme  to  update  all  cells, 

«r'  =u?-^-(fi+u2-f.v2) 

at  the  regular  cells,  i* 0,  and 

u0+1  ~  “0  “  ( f\l2  -f-lll) 

at  the  small  cell.  We  want  to  define  the  fluxes /±1/2  so  that  the  overall  scheme  is  stable  as  a  -»0. 

Typical  flux  functions  for  the  regular  cells  include  Godunov’s  method,  which  for  a  scalar  equations 
with/'(u)  5  0  is  just  upwind  differencing,  with 

fi+V2  -F(ui’Un.\)  =/(«i+l) 

and  the  second  order  Lax  Wendroff  scheme 

(U,+Uj+ 1)  A  t 

F(u„Ui+l)=  -2 —  +  *,)-/(«,)). 

However,  these  definitions  must  be  modified  at  the  small  cell,  since  if  we  define  fv2  =  F(u0,ut),  the 
resulting  scheme  loses  accuracy  and  becomes  unstable  for  small  a.  Our  approach  is  to  define  a  new  state 
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Figure  1  Special  flux  formulas  arc  needed  at  the  edges  of  the  small  cell  to  maintain  stability  and  accuracy. 


uitf,  and  let 


fv2=F(u,,f„Ul). 

Here,  u,,f,  is  an  approximation  to  u  at  distance  hi 2  from  the  interface.  For  example,  uuf,  might  be  obtained 
by  linear  interpolation  between  u_,  and  u0, 

2au0  +  (l-a)u_! 


Note  that  ultf,~*u.0  as  a-»l,  and  uUfl—>u  j  as  a-»0.  In  either  of  these  limits  the  grid  becomes  regular  again 
and  the  difference  scheme  reverts  to  the  uniform  scheme. 

Similarly,  we  define  a  flux/_1/2  using  the  left  state  u. 1  and  a  new  right  state  un(hl  defined  by  interpo 
lating  u  to  a  point  a  distance  hl2  to  the  right  of  this  interface,  e.g. 

2au0  +  (l-a)ui 

bright  ~  "I  ■ 

*  1  +  a 

We  then  set 


f-V2  =  ^(«-l 

It  can  be  shown  for  the  equation  u,  =  u,  that  if  uttf,  and  un/ht  are  obtained  by  linear  interpolation  then 
the  resulting  scheme  is  stable  as  a-»0,  using  a  time  step  Ar  that  satisfies  the  CFL  condition  for  the  regular 

grid,  This  results  holds  for  both  upwind  differencing  and  Lax  WendrofF.  However,  for  Lax  Wendr- 
n 

off,  a  more  accurate  procedure  would  be  to  use  quadratic  interpolation  for  ultf,  and  ung>u.  In  this  case,  a 
combination  of  theoretical  and  numerical  results  show  that  if  the  additional  point  used  in  the  interpolation 
is  the  upwind  point  then  the  scheme  is  stable  as  a— >0,  but  if  the  downwind  point  u_j  is  used  ,  then  the  CFL 
limit  is  reduced  to  1/2. 


Stable  Difference  Equations  for  Two  Dimensional  Irregular  Cells 

For  two  dimensional  calculations,  we  again  need  a  new  difference  scheme  to  compute  fluxes  at  the 
edges  of  the  irregular  cells  adjacent  to  the  boundary.  The  boundary  of  a  solid  body  is  represented  by  a 
piecewise  linear  segment  in  each  cell,  so  that  the  irregular  cells  can  have  either  3,  4  or  5  sides.  Our 
approach  uses  locally  normal  and  tangential  coordinate  directions  to  define  left  and  right  states  for  a 
Riemann  problem  at  each  cell  edge.  This  fits  naturally  with  the  MUSCL  scheme  used  in  the  interior  of  the 
domain  in  our  calculations.  However,  we  generalize  the  pointwise  approach  in  the  one  dimensional  case, 
and  use  aconservative  averages  of  the  solution  over  a  box  a  distance  h  in  the  appropriate  direction  away 
from  a  cell  edge..  For  example,  in  Figure  2  the  state  qt  is  obtained  using  an  area-weighted  average  of  the 
values  Uij-i  and  j  that  intersect  the  box  from  the  regular  grid.  In  an  analogous  way  we  obtain  the  state 
qr.  These  values  are  then  rotated  into  a  frame  of  reference  that  is  tangent  to  the  boundary,  and  a  one  dimen¬ 
sional  Riemann  problem  in  the  tangential  direction  is  solved.  This  gives  the  flux  /r.  This  procedure  is 
repeated  for  the  dashed  boxes  in  Figure  2  that  are  normal  to  the  boundary,  giving  a  flux  fA  in  the  t)  direc¬ 
tion.  (The  part  of  the  box  that  lies  outside  the  domain  is  interpolated  from  qk,  see  Figure  3).  The  final  value 
of  the  flux  at  the  vertical  interface  is  a  linear  combination  /$ cosG  +  /nsin0,  where  0  is  the  angle  the  boun¬ 
dary  makes  with  the  grid.  For  a  boundary  with  curvature,  we  determine  these  directions  using  the  boun¬ 
dary  segment  of  the  cell  with  the  smaller  area  adjacent  to  the  interface.  This  helps  retain  stability  for  the 


Figure  2  shows  a  schematic  of  the  rotated  difference  scheme  used  to  define  the  vertical  flux. 


smaller  cells,  by  maintaining  a  certain  cancellation  property  of  our  flux  definitions,  described  more  com¬ 
pletely  in  [Berger  and  LeVeque].  Related  work  using  rotated  difference  schemes  has  been  done  by  [Jame¬ 
son;  Davis;  Levy,  Powell  and  van  Leer], 

At  the  solid  wall  boundary  itself,  the  flux  can  be  determined  more  simply,  using  only  boxes  normal 
to  the  boundary  as  shown  in  Figure  l.  First  we  obtain  a  value  qk  for  the  box  interior  to  the  domain,  using 
area  weighted  averages,  and  rotate  the  velocities  into  the  boundary  coordinate  frame.  A  boundary  Riemann 
problem  is  solved  between  qk  and  qk  (with  negated  normal  velocity),  to  satisfy  the  the  boundary  conditions 
of  no  normal  flow. 


Figure  3  indicates  the  scheme  used  to  determine  the  boundary  flux. 

In  this  case,  if  the  solid  wall  boundary  happens  to  align  with  the  Cartesian  grid,  the  scheme  reverts  to 
the  usual  first  order  Godunov  method.  To  improve  the  scheme  to  second  order,  following  the  MUSCL 
approach  as  described  in  [Colella],  we  need  to  introduce  limited  slopes  in  the  solution  reconstruction  phase, 
and  tangential  derivatives  for  predicting  states  at  the  cell  edges.  These  steps  are  also  necessary  to  improve 
the  stability  limit  for  Godunov’s  method  from  1/2  to  1.  Work  on  these  improvements  is  continuing.  Refer¬ 
ring  to  Figure  3,  we  add  a  tangential  derivative  /$  to  the  state  qk  for  the  normal  Riemann  problem,  with 
f\-f  (u(Jc,k+l))-f  (u(k-l,k))  /  /*,  where  /*  is  the  length  of  the  k*  boundary  segment.  The  state 
u(k,k  + 1)  comes  from  solving  a  Riemann  problem  in  the  tangential  direction  at  the  interface  between  cells  k 
and  k+ 1.  As  before,  the  stencil  for  this  Riemann  problem  must  be  enlarged  beyond  the  adjacent  cells  to 
maintain  stability.  For  example,  the  right  state  at  this  interface  is  not  just  the  value  qk +I ,  but  a  linear 


* 


combination  of  the  solution  in  several  boxes,  qk+l  and  qk*2<  up  to  a  distance  h  away  from  the  interface. 
The  left  state  at  this  right  interface  can  be  taken  to  be  the  value  qk  since  the  length  of  that  cell's  boundary 
segment  is  larger  than  h.  This  same  procedure  is  used  to  include  tangential  derivatives  in  the  normal  box 
Ricmann  problems  for  interior  cell  edges.  This  procedure  alone  improves  the  CFL  limit  to  1.  It  remains  to 
incorporate  monotonized  slopes  into  the  scheme  in  order  to  achieve  second  order  accuracy. 

While  the  overall  scheme  at  the  boundary  involves  twice  as  many  Riemann  problems  as  the  ordinary 
MUSCL  scheme,  it  is  fully  vectorizable.  The  coefficients  in  the  interpolations  for  the  left  and  nght  states 
are  fixed  for  the  duration  of  the  integration,  and  are  not  dependent  of  the  properties  of  the  solution  at  each 
step.  In  numerical  experiments  in  two  dimensions,  this  scheme  remains  stable  for  cell  areas  that  are  orders 
of  magnitude  smaller  than  the  regular  cell  areas  (down  to  the  round-off  level).  In  essence,  our  method  can 
be  viewed  as  a  technique  for  defining  fluxes  on  an  irregular  grid  by  a  very  local  mapping  to  a  regular  grid. 
This  viewpoint  may  prove  useful  in  defining  higher  order  methods  on  unstructured  grids. 

Computational  Example 

We  illustrate  this  by  computing  time  dependent  flow  around  a  cylinder.  The  initial  conditions  are  an 
incident  shock  traveling  at  Mach  2.81.  We  use  a  simple  MUSCL  scheme  to  advance  the  flow  field  in  the 
interior  of  the  domain.  Figure  4  shows  a  contour  plot  of  the  flow  field,  as  well  as  a  plot  of  density  as  a  func¬ 
tion  of  arclength  around  the  cylinder.  The  only  cells  drawn  on  the  contour  plot  are  the  irregular  cells  from 
the  Cartesian  grid  that  intersect  the  body.  Note  the  smoothness  of  the  arclength  plot  despite  the  irregularity 
of  the  grid  around  the  body.  This  example  was  computed  using  the  local  mesh  refinement  of  [Berger  and 
Colella];  the  location  of  the  rectangular  fine  grids  is  indicated  on  the  contour  plot  as  well. 
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ABSTRACT 

We  describe  a  new  point  clustering  algorithm,  and  its  application  to  automatic  grid 
generation,  a  technique  used  to  solve  partial  differential  equations.  Algorithms  from  the 
computer  vision  and  pattern  recognition  literature  are  used  to  partition  points  into  a  set  of 
enclosing  rectangles.  We  show  examples  from  two  dimensional  calculations,  but  the 
algorithm  generalizes  readily  to  three  dimensions. 


1.  Introduction 

This  paper  presents  a  new  point  clustering  algorithm  and  its  application  to  automatic  grid  generation. 
The  algorithm  clusters  the  points  into  distinct  rectangles  such  that  neighboring  points  are  in  the  same  rec¬ 
tangle  (as  much  as  possible),  and  all  points  are  contained  in  some  rectangle.  The  application  is  best  illus¬ 
trated  by  an  example.  We  are  solving  a  partial  differential  equation  using  finite  difference  techniques.  The 
difference  equations  are  first  solved  on  the  uniform  coarse  grid  in  Figure  1.  An  "error  estimation"  pro¬ 
cedure  [Oliger]  is  then  used  to  flag  grid  points  that  need  to  be  in  a  finer  grid,  which  is  just  a  smaller  rec¬ 
tangular  grid  with  finer  mesh  spacing.  Figure  1  shows  the  flagged  coarse  grid  points  as  well  as  the  new  fine 
grids  that  together  contain  all  the  flagged  points.  This  procedure  is  usually  employed  recursively,  leading 
to  a  nested  sequence  of  increasingly  fine,  locally  uniform  subgrids  which  are  superimposed  on  the  underly¬ 
ing  base  grid. 

The  grid  generation  problem  then  is  to  define  a  set  of  rectangles  that  enclose  all  the  flagged  points. 
Certain  factors  make  a  huge  difference  in  the  performance  of  these  grids  in  solving  the  pdes. 

[1]  There  should  be  as  Utile  unnecessarily  refined  area  as  possible. 

Since  the  cost  of  the  numerical  integration  procedure  is  proportional  to  the  area  of  the  rectangle,  the 
smaller  the  better.  Figure  1  gives  an  example  of  a  set  of  flagged  points  for  which  a  few  patches  lead  to 
much  less  refined  area  than  refining  the  whole  grid.  This  gain  in  efficiency  is  the  purpose  of  adaptive 
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Figure  1  shows  a  coarse  grid  with  rectangular  subgrids  around  the  flagged  points.  The  coarse  grid 
points  with  high  error  are  marked  with  an  "x". 

methods.  Some  unnecessarily  refined  area  (or  inclusion  of  non-fiagged  coarse  grid  points  in  a  new  rectan¬ 
gle)  is  inevitable,  since  we  are  restricted  to  using  rectangles.  In  addition,  for  numerical  reasons  the  rectan¬ 
gles  are  oriented  with  the  base  grid  rectangle.  This  is  true  even  if  the  flagged  points  lie  on  a  diagonal  of  the 
coarse  grid,  and  could  be  perfectly  enclosed  by  a  rotated  rectangle.  (However,  an  algorithm  that  uses 
rotated  rectangles  is  considered  in  [Berger]).  Along  these  lines,  if  several  rectangles  are  used  to  enclose  the 
flagged  points,  their  overlap  should  be  minimal. 

Another  criterion  for  generating  these  rectangles  is: 

[2]  There  should  be  as  few  rectangles  as  possible. 

At  the  other  extreme,  we  could  put  one  tiny  rectangle  around  each  flagged  point.  Many  of  these  tiny  rec¬ 
tangles  would  share  a  common  boundary  segment  However,  there  is  boundary  overhead  associated  with 
each  rectangular  subgrid  that  should  also  be  minimized,  along  with  the  area.  In  addition,  these  procedures 
will  be  used  on  vector  processors,  which  favors  larger  vector  lengths  and  therefore  larger  rectangles.  (We 
could  worcy  further  about  this,  for  example  by  trying  to  maximize  the  length  in  a  particular  coordinate 
direction,  but  we  will  not  consider  such  machine  specific  details  here). 
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The  third  criterion  is  the  most  difficult  to  pinpoint 

[3]  The  rectangles  should  "fit"  the  data. 

This  is  hard  to  make  absolutely  precise,  but  for  example,  if  a  person  were  to  put  the  rectangles  around  the 
points  by  hand,  using  whatever  clustering  or  partition  the  brain  uses,  it  would  "look  right".  Although  this  is 
not  essential  for  accurate  numerical  integration  on  the  rectangles,  we  prefer  the  adaptively  generated 
subgrids  to  be  pleasing.  Finally, 

[4]  The  algorithm  should  be  fast. 

This  algorithm  is  repeated  every  few  timesteps,  or  hundreds  of  times  during  any  particular  numerical  simu¬ 
lation,  and  should  therefore  be  fast  relative  to  the  time  needed  to  take  a  time  step  on  the  resulting  grids. 

Our  solution  to  this  rectangle-fitting  problem  uses  algorithms  from  the  pattern  recognition  and  com¬ 
puter  vision  literature.  A  combination  of  signature  arrays  and  zero  crossings  of  second  derivatives  is  used 
to  partitnn  the  flagged  points  into  rectangles.  Our  examples  are  all  in  two  dimensions,  but  the  algorithm 
generalizes  readily  and  has  proven  effective  in  three  dimensions  too.  Before  describing  our  algorithm,  we 
give  a  little  background  and  discuss  some  other  approaches  we  tried  and  discarded. 

1.1.  Previous  Algorithms 

Our  previous  algorithms  for  this  problem  can  be  summarized  as  being  of  two  main  types:  bottom-up 
or  top-down.  The  top-down  approach  is  based  on  a  bisection  method,  it  can  be  viewed  as  a  form  of 
divisive  hierarchical  clustering  [Duda  and  Hart].  Initially,  the  flagged  points  of  the  grid  are  surrounded  by 
a  single  rectangle,  and  its  efficiency  is  computed.  Here,  we  define  the  efficiency  of  a  grid  as  the  ratio  of 
flagged  points  to  the  total  number  of  coarse  grid  points  in  the  new  rectangle.  This  is  one  of  the  key  parame¬ 
ters  behind  our  algorithm,  and  it  is  easily  computed.  If  the  efficiency  is  above  a  preselected  threshold,  the 
rectangle  is  accepted  and  the  algorithm  stops.  Otherwise,  we  bisea  the  longest  direction  of  the  rectangle, 
and  generate  two  new  subgrids.  This  process  is  repeated  recursively  on  each  of  the  two  subgrids.  When  the 
algorithm  terminates,  all  of  the  subgrids  are  guaranteed  to  have  acceptable  efficiency  ratings.  However, 
hierarchical  clustering  methods  are  known  to  create  clusters  even  if  no  natural  clusters  exist  [Anderberg; 
Haitigan;  Jain  and  Dubes].  In  addition,  since  the  bisection  uses  no  information  about  the  locations  of  the 
flagged  points,  a  non -optimal  grid  hierarchy  is  generally  created.  To  alleviate  this,  we  usually  follow  the 
bisection  step  with  a  merging  step,  where  neighboring  subgrids  are  merged  into  larger  subgrids  if  the  result 
continues  to  be  acceptably  efficient  This  mrrging  step  is  what  to  the  problem  of  overlapping  grids. 
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The  bottom -up  approach  starts  at  the  grid  point  leveL  The  flagged  points  are  organized  into  a 
minimal  spanning  tree,  so  that  each  point  is  connected  to  its  nearest  neighbor.  Neighboring  branches  of  the 
tree  are  merged,  either  one  point  at  a  time,  or  a  front  at  a  time,  as  long  as  the  resulting  grid  is  efficient 
Although  philosophically  appealing,  this  algorithm  actually  performs  much  worse  than  the  top-down  bisec¬ 
tion  algorithm.  A  fundamental  problem  is  the  non-uniqueness  of  the  minimal  spanning  tree.  Also,  the 
algorithm  suffered  from  the  hill-climbing  problem  of  getting  stuck  in  local  minima;  it  was  very  sensitive  in 
the  beginning  steps  of  the  algorithm  to  the  initial  direction  of  growth  of  the  clusters,  and  tended  to  stop 
prematurely,  although  larger  and  acceptably  efficient  grids  were  just  several  branches  away.  In  that  case,  it 
had  to  be  followed  by  a  merging  procedure  as  well. 
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Figure  2  Nearest  neighbor  clustering  is  sufficient  in  simple  cases,  e.g.,  the  top  left  cluster. 


In  practice,  both  approaches  were  preceded  by  a  "nearest  neighbor"  clustering  algorithm.  The  pur¬ 
pose  of  this  was  to  separate  flagged  points  when  possible  into  isolated  islands  (see  Figure  2).  This  some¬ 
times  produces  acceptable  clusters  by  itself,  but  fails  to  help  when  the  flagged  points  formed  elongated, 
curved  shapes.  Thus,  it  was  followed  by  either  the  bisection  or  minimal  spanning  tree  algorithm.  Summar¬ 
izing,  both  of  these  algorithms  produced  less  than  optimally  efficient  grids  that  overlapped  too  much. 
Better  grids  were  easily  created  by  hand. 

2.  Towards  an  Efficient  Algorithm 

In  a  more  general  form,  die  grid  generation  algorithm  should  cluster  a  set  of  m  flagged  points  into  k 
clusters,  where  k  is  either  specified  a  priori  or  is  determined  by  the  algorithm  itself.  This  special  type  of 
clustering  is  called  partidoning  [Anderberg],  Our  first  approach  considered  the  question  of  how  to  choose  a 
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set  of  k  "seed  points"  around  which  the  k  clusters  would  be  built 

2.1  A  First  Approach:  Local  Maxima  in  Two-Dimensional  Grids 

Initially,  each  of  the  grid  points  is  given  a  value:  "1”  for  the  flagged,  and  "0"  for  the  non-flagged 
ones.  These  grid  values,  viewed  as  a  binary  image,  /(x,y),  are  then  preprocessed  by  convolving  it  with 
either  an  "averaging"  or  a  "low-pass"  filtering  template,  (see  Figure  3)  [Ballard  and  Brown;  Horn;  Levine], 

_  A 

This  operation  results  in  a  non-convex  function,  /(x,y),  whose  local  maxima,  determined  by  the  Sobel 
operators  [Ballard  and  Brown;  Levine]  of  Figure  4,  compose  the  set  of  "seed  points”  around  which  we 
build  the  clusters. 


b) 


Figure  3  shows  the  filtering  templates:  (a)  Averaging  (b)  Low-Pass. 

Three  partitioning  algorithms  were  tested  using  the  seeds  found  above:  the  standard  k-means  algo¬ 
rithm,  its  converging  variant,  and  a  k-means  variant  where  no  updating  of  the  centroids  takes  place 
[Anderberg;  Hartigan].  These  algorithms  are  outlined  in  the  Appendix. 
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Figure  4  shows  the  Sobel  gradient  operators:  (a)  d/dx  (b)  3/d y. 


Figures  S  and  6  show  graphically  the  output  of  the  three  algorithms  on  two  sample  data  sets.  The 
three  algorithms  exhibit  the  same  overall  behavior.  The  seeds  are  sensibly  chosen,  but  in  all  the  test  cases 
the  resulting  grids  overlap  excessively.  We  attribute  this  problem  to  the  large  number  of  seeds  discovered 
by  this  method.  The  next  method  tries  to  reduce  the  number  of  seeds,  keeping  only  the  best,  and  thus  hope¬ 
fully  reducing  the  overlapping. 


Figure  5  The  seeds  are  local  maxima  of  a  two  dimensional  function. 


J 


Figure  6  The  seeds  are  local  maxima  of  a  two  dimensional  function. 
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22  A  Second  Approach:  Local  Maxima  in  Signature  Arrays 

Signatures  have  been  known  in  computer  vision  and  pattern  recognition  for  many  years  [Ballard  and 
Brown;  Horn;  Levine].  Able  to  encapsulate  gross  characteristics  of  a  shape,  and  computationally  simple, 
signatures  proved  very  useful  for  establishing  preliminary  landmarks  in  images;  these  landmarks  in  turn  led 
to  a  subsequent  reduction  of  the  search  effort.  Given  a  continuous  function,  the  horizontal  and  vertical  sig¬ 
natures,  H(x)  and  V(y)  are  defined  as 

W(x)  =  J/(x,y)dy 

y 

and 

V(y)  =  \f(x,y)dx 

X 

respectively. 

First,  the  horizontal  and  vertical  signatures  of  the  image  are  computed.  The  resulting  one¬ 
dimensional  arrays  are  "smoothed"  [Horn]  using  the  template  in  Figure  7,  and  subsequently  searched  for 
local  maxima.  Let  M  and  N  be  the  two  sets  of  maxima.  After  discarding  those  tuples  of  the  Cartesian  pro¬ 
duct  MxN  that  correspond  to  non-flagged  regions,  we  are  left  with  precisely  the  coordinates  of  the  starting 
seeds. 


1  1  2  |  1 

Figure  7  shows  the  one-dimensional  smoothing  template. 

With  this  choice  of  seeds,  we  again  employed  the  three  partitioning  techniques  (k-means,  converging 
k-means,  and  the  no  updating  variant).  Figures  8  and  9  show  the  results.  This  algorithm  resulted  in  fewer 

A 

subgrids  and  reduced  overlapping,  making  this  approach  superior  to  using  the  local  maxima  of  I(x,y)  for 
the  seed  points.  Unfortunately,  the  generated  subgrids  were  still  not  the  most  efficient  ones;  better  choices 
were  clearly  possible.  Some  observations  based  on  extensive  experimentation  could  still  be  made:  1)  the 
non-converging  variants  outperformed  the  converging  k-means,  2)  overlapping  was  minimal  when  no 
updating  of  the  centroids  occurred,  and  3)  the  distance  metric  (Manhattan  Block  vs.  Euclidean)  had  no 
appreciable  effect  on  the  results. 

Apparently,  the  use  of  local  maxima  in  signatures  did  not  capture  enough  of  the  underlying  structure. 
In  the  next  section,  we  use  signatures  in  a  different  way  to  partition  the  flagged  points  into  clusters. 

3.  The  Algorithm 

Our  best  algorithm  uses  ideas  related  to  edge  detection.  One  of  the  many  approaches  to  the  edge 
detection  problem  is  the  one  suggested  by  [Marr  and  Hildreth].  Based  on  the  psychophysical  and  neuro¬ 
physiological  experiments  of  (Campbell  and  Robson],  the  method  consists  of  first  convolving  the  original 


Figure  9  The  seeds  are  chosen  from  local  maxima  of  signature  arrays. 
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image  against  a  Gaussian  kernel  and  then  computing  the  Laplacian  of  the  result;  the  intensity  discontinui¬ 
ties  are  associated  with  those  positions  where  the  Laplacian  is  equal  to  0  (zero  crossings). 

In  what  follows,  the  input  grids  are  viewed  as  binary  images  in  the  sense  of  paragraph  2.1.  The 
edges  will  now  be  located  at  those  positions  of  the  grid  where  a  transition  from  a  flagged  point  region  to  a 
non-flagged  one  occurs.  The  most  prominent  such  transition  represents  a  “natural"  line  with  respect  to 
which  the  original  grid  can  be  partitioned.  For  the  example  depicted  in  Figure  10,  the  line  (e)  represents 
such  a  transition. 


X  X 

(e)  X  X 

X  X  X  X 

X  X  X  X 
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X  X  X  X  X 

Figure  10  The  rectangle  is  partitioned  at  line  (e);  after  this,  an  isolated  cluster  can  be  detected  for 
further  partitioning. 


The  problem  is  that  of  determining  such  a  transition.  The  Marr-Hildreth  method,  although  a  poten¬ 
tial  candidate,  cannot  easily  be  employed  towards  this  end  for  two  reasons.  First,  its  output  is  a  collection 
of  Boolean  values  at  the  different  grid  locations:  TRUE  if  a  zero  crossing  exists  at  that  location,  FALSE 
otherwise;  these  Boolean  values  need  to  be  combined  in  order  to  generate  an  actual  edge,  a  not  so  trivial 
task.  Second,  since  the  Laplacian  operator  is  isotropic,  the  “natural"  line  with  respect  to  which  the  original 
grid  could  be  split  will  not,  in  general,  be  parallel  to  the  sides  of  the  grid. 

Signatures  again  hold  the  answer  the  idea  is  to  look  for  zero  crossings  in  the  second  derivative  of  a 
signature.  This  idea  borrows  from  both  the  signature  approach  and  the  Marr-Hildreth  method,  except  that 
we  do  not  convolve  with  a  Gaussian  filter. 

In  general,  there  is  more  than  one  zero  crossing  in  a  given  signature  array.  For  our  purposes  how¬ 
ever,  we  only  select  one  zero  crossing  at  a  time;  it  corresponds  to  most  prominent  edge,  and  its  location  is 
determined  by  searching  both  the  horizontal  and  vertical  signature  arrays  for  the  zero  crossing  correspond¬ 
ing  to  the  largest  local  change  of  values.  This  is  indicated  in  Figure  11,  where  the  row  (resp.  column) 
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Figure  11  The  steepest  zero  crossing  occurs  at  line  (e),  and  makes  the  most  efficient  rectangles. 

labelled  Z  contains  the  horizontal  (vertical)  signature,  and  the  row  (column)  labelled  A  contains  the  Lapla- 
cian  of  the  signature  in  the  appropriate  direction. 

The  input  grids  are  also  expected  to  contain  isolated  regions  of  flagged  points  (islands),  making  it 
necessary  that  both  signature  arrays  be  first  searched  for  chains  of  0’s,  or  “holes".  This  can  be  seen  in  Fig¬ 
ure  10.  The  occurrence  of  such  holes  provides  obvious  choices  for  splitting  the  input  grid  into  a  number  of 
subgrids,  and  is  exploited  before  any  attempt  at  locating  zero  crossings  is  made.  If  any  holes  are  found,  we 
do  not  attempt  to  compute  any  zero  crossings.  We  proceed  by  applying  the  above  steps  to  only  those  of  the 
generated  subgrids  that  are  still  inefficient. 

We  now  give  a  high  level  description  of  the  actual  algorithm,  followed  by  sample  runs  in  Figures  12 
to  16  on  a  number  of  real  and  synthetic  problem  cases. 
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* 


BEGIN 

i-1; 

while  ( i  <=  number_of_rec tangles  )  do 

if  (rectangle  efficiency  <  threshold )  then 
compute  signatures ; 

find  the  best  place  to  split  ( either  a  hole  or  edge ) ; 
if  ( found  a  place  to  split )  then 
split  rectangle  in  two  ; 

append  new  rectangle  to  end  of  list  of  rectangles  ; 
else 

i  =  i+  1 ; 
endif 
else 

i  =  i  +  1 ;  ( consider  the  next  rectangle  on  the  list ) 
endif 
end  while 

END 


4.  Comments  on  the  Algorithm  and  its  Performance 

While  most  of  the  time  the  algorithm  performs  exceptionally  well,  sometimes  it  generates  anomalous 
and/or  non -optimal  rectangles.  As  we  will  see,  many  of  these  can  be  eliminated  by  simple  changes  in  the 
algc  ruhm.  The  anomalies  fall  into  several  general  categories,  which  we  illustrate  by  picture.  The  basic 
algorithm  description  has  two  exit  points:  a  rectangle  is  "accepted"  either  because  its  efficiency  is  already 
above  threshold,  or  because  it  cannot  be  further  split  using  either  of  the  two  methods  (Le.  holes  and 
inflection  points).  As  a  result,  the  efficiency  of  some  of  the  generated  subgrids  may  be  below  the  preset 
threshold.  This  is  the  case  with  all  grids  associated  with  regions  that  form  an  angle  with  the  horizontal.  The 
inefficiency  approaches  its  maximum  as  angles  approach  45  degrees. 

Another  problematic  grid  is  the  one  appearing  in  Figure  17:  as  can  be  seen  neither  the  horizontal  nor 
the  vertical  signatures  will  contain  any  holes  or  zero  crossings,  and  this  will  be  true  for  all  non-zero  values 
of  ajo.  In  all  such  cases  the  bounding  rectangles  will  have  an  efficiency  of  precisely  50%.  An  ordinary 
bisection  step  could  be  easily  incorporated  here  to  increase  the  efficiency  to  100%.  A  similar  arrangement 
of  flagged  points  which  we  have  encountered  in  our  experiments  (see  Figure  18)  leads  to  another  kind  of 
non-optimal  choice.  One  way  around  this  is  to  used  weighted  second  derivatives,  scaling  the  Laplacian  by 
the  number  of  flagged  points.  This  leads  to  a  correct  choice  for  the  zero  crossing  in  Figure  18. 


Figure  13  An  example  of  die  final  algorithm. 


Figure  15  An  example  of  the  final  algorithm. 


r  xxxx’m'x^i 

;  xxxxxxxxx 

;  xxxxxxxxx 

xxxxxxxxx  b 
xxxxxxxxx! 

!  ^ _ f _ ^  XXXXXXXXX 

xxxxxxxxx  ** - a - •*  ; 

XXXXXXXXX  ! 

xxxxxxxxx 

4  'XXXXXXXXX  ! 

Xxxxxxxxx  ! 

tomxxxx _ j 


Figure  17  There  are  no  zero  crossings  in  the  signature  arrays. 
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Figure  18  This  set  of  points  leads  to  a  non-optimal  partition,  unless  the  zero  crossings  are  scaled. 


A  modification  of  our  algorithm  that  covers  both  anomalous  cases  is  to  compute  the  sum  of  the  abso¬ 
lute  value  of  the  gradient,  and  difference  the  results  to  get  the  second  derivative.  The  most  robust  solution 
to  this  problem  is  still  an  open  question. 

A  seemingly  problematic  case  is  shown  in  Figure  19.  Here  it  appears  as  if  the  algorithm  made  a 
non-optimal  decision  by  unnecessarily  splitting  rectangle  R 1  (dotted  line)  into  two  smaller  subrectangles. 


! 
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However,  careful  inspection  shows  that  rectangles  R 1  and  R  2  could  not  have  been  generated  without  intro¬ 
ducing  an  overlapping  region. 


Figure  19  shows  an  unexpected  set  of  rectangles. 

Tight  bounds  for  the  running  time  of  the  described  algorithm  are  very  hard  to  establish,  since  the  pre¬ 
cise  flow  of  the  algorithm  is  input  dependent.  However,  it  should  be  clear  that  the  running  time  is 
0(k(P+N+M)),  where  A:  is  the  total  number  of  grids  upon  termination  of  the  algorithm,  and  P  is  the 
number  of  flagged  points.  The  P  term  comes  from  computing  the  signatures  of  the  points  (by  traversing  a 
list  of  the  flagged  points).  The  N  (resp.  Af)  tom  comes  from  the  linear  search  that  determines  the  best 
inflection  point  The  above  bound  is  by  no  means  optimal,  and  the  algorithm  performs  very  well  in  prac¬ 
tice.  Preliminary  three  dimensional  results  also  show  great  performance. 

5.  Conclusion 

We  have  described  a  new  and  efficient  algorithm  for  point  clustering  and  adaptive  grid  generation. 
Hie  algorithm's  performance  has  been  demonstrated  through  a  series  of  graphs  showing  results  obtained 
with  both  synthetic  and  actual  2-dimensional  inputs.  Preliminary  experiments  with  3-dimensional  problems 
also  show  a  considerably  improved  performance  over  the  previous  approaches.  In  general,  the  efficiency  of 
the  enclosing  rectangles  for  our  applications  has  been  very  high,  typically  ranging  between  83%  and  100%, 
with  the  exception  of  the  problematic  cases  illustrated  in  figures  17-19.  It  is  surprising  how  effective  the 
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algorithm  performs  on  multidimensional  data,  even  though  it  is  based  on  Cartesian  coordinate  directions. 

Finally,  this  algorithm  may  also  prove  useful  in  other  applications  with  binary  image  data,  for  example  in 

generating  bounding  rectangles  for  computer  graphics  applications.  A  rectangle  fitting  algorithm  has  also 

been  used  in  conjunction  with  a  pattern  recognition  system  for  understanding  Japanese  business  cards  [Kise 

et  al.]. 

6.  Appendix 

6.1  McQueen’s  1-means  Partitioning  Algorithm 

[1]  Form  1  single  member  clusters  each  one  containing  precisely  one  of  the  k  starting  seeds.  The  clus¬ 
ters’  centroids  originally  coincide  with  the  starting  seeds. 

[2]  Assign  each  of  the  remaining  data  points  to  the  cluster  with  the  nearest  (with  respect  to  an  appropri¬ 
ate  distance  metric)  centroid  recomputing  the  gaining  cluster’s  centroid  after  each  assignment. 

[3]  Assume  the  cluster  centroids  are  fixed  this  time,  and  reassign  each  of  the  data  points  to  the  cluster 
with  the  nearest  centroid  (one  pass  through  the  data). 

(J  McQueen’s  Converging  1-means  Partitioning  Algorithm 

[1]  FOrm  k  single  member  dusters  each  one  containing  precisely  one  of  the  1  starting  seeds.  The  clus¬ 
ters’  centroids  originally  coincide  with  the  starting  seeds. 

[2]  Assign  each  of  the  remaining  data  points  to  the  duster  with  the  nearest  (with  respect  to  an  appropri¬ 
ate  distance  metric)  centroid  recomputing  the  gaining  cluster’s  centroid  after  each  assignment. 

[3]  For  each  data  point  compute  its  distance  to  all  the  cluster  centroids;  if  the  nearest  centroid 
corresponds  to  a  cluster  other  than  the  point’s  actual  parent  duster  reassign  the  point;  recompute  the 
centroids  of  both  the  gaining  and  losing  clusters. 

[4]  Repeat  step  3)  until  a  full  sweep  through  the  data  does  not  induce  further  changes  in  the  points’ 
memberships. 


63  Non-Updating  Variant  of  1-means  Partitioning  Algorithm 

[1]  Form  1  single  member  dusters  each  one  containing  precisely  one  of  the  1  starting  seeds.  The  clus¬ 
ters’  centroids  remain  fixed  throughout  the  algorithm. 

[2]  Assign  each  of  the  remaining  data  points  to  the  duster  with  the  nearest  (with  respect  to  an  appropri¬ 
ate  distance  metric)  centroid  but  do  not  update  the  gaining  cluster’s  centroid  (one  pass  though  the 
data). 
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Abstract— The  inviscid  Euler  equations  in  complicated  geometries  are  solved  using  a  Cartesian  grid.  This 
requires  wall  boundary  conditions  in  the  irregular  grid  cells  near  the  boundary.  Since  these  cells  may  be 
orders  of  magnitude  smaller  than  the  regular  grid  cells,  stability  is  a  primary  concern.  A  new  approach 
to  this  problem  is  presented  and  illustrated. 


t.  INTRODUCTION 

In  previous  work1  i  we  have  described  a  Cartesian 
grid  method  for  the  inviscid  Euler  equations  in 
arbitrary  gemetries.  There  are  many  advantages  to 
be  gained  from  this  approach.  Grid  generation  is 
simplified,  since  we  avoid  the  use  of  (possibly  multi¬ 
block)  body-fitted  grids,  and  we  can  use  high  resol¬ 
ution,  highly  efficient  solvers  on  regular  grids  over  the 
bulk  of  the  domain.  This  has  led  to  renewed  interest 
in  Cartesian  grids  in  recent  years  (see  Refs  4  and  5). 
One  of  the  difficulties  with  Cartesian  grids  is  that 
they  give  insufficient  resolution  in  certain  regions 
such  as  leading  edges.  This  can  now  be  overcome  by 
Cartesian  adaptive  mesh  refinement.1-6 

The  principal  remaining  difficulty  in  this  approach 
is  due  to  the  essentially  arbitrary  way  that  a  Cartesian 
grid  intersects  the  boundaries  of  the  computational 
domain.  In  particular,  a  solid  wall  boundary  cutting 
through  the  grid  creates  irregular  cells  that  may  be 
orders  of  magnitude  smaller  than  the  regular  cells 
away  from  the  boundary.  For  these  irregular  cells, 
special  difference  equations  are  needed  that  maintain 
stability  and  accuracy,  and  satisfy  the  solid  wall 
boundary  conditions  of  no  normal  flow. 

In  this  work,  we  present  an  improved  method  for 
the  small  boundary  cells.  We  use  an  explicit,  finite 
volume  formulation  that  computes  fluxes  at  cell  edges 
on  the  regular  part  of  the  domain.  We  would  like  to 
define  fluxes  at  the  edges  of  the  irregular  cells  in  such 
a  way  that  the  method  is  stable  even  with  very  small 
boundary  cells,  using  a  time  step  based  on  the  regular 
grid  cells  away  from  the  boundary.  The  Courant- 
Friedrichs-Lewy  (CFL)  condition  requires  that  the 
numerical  method  allows  information  to  propagate  at 
least  as  quickly  as  the  underlying  differential  equation. 
In  the  present  context  this  means  that  we  must  define 
fluxes  at  the  sides  of  our  irregular  cells  based  on  more 
than  just  the  neighboring  cell  values. 

In  our  previous  work,  we  have  used  a  wave 
propagation  approach  in  defining  these  fluxes  Here 


we  propose  an  alternative  method  that  has  some 
advantages  over  the  wave  propagation  approach.  In 
particular,  the  wave  propagation  method  is  subject 
to  intermittent  instabilities  due  to  two-dimensional 
effects  that  are  not  clearly  understood.  The  new 
method  has  a  cancellation  property  in  two  dimensions 
that  appears  to  give  better  stability  properties.  More¬ 
over,  the  computational  geometry  is  simplied  in  the 
new  approach.  The  fluxes  are  defined  in  terms  of 
weighted  averages  of  nearby  cell  values.  These  weights 
may  be  calculated  as  a  preprocessing  step  on  any 
fixed  grid  and  need  not  be- repeatedly  calculated.  In 
the  previous  approach  the  weights  depended  on  the 
flow  variables  and  a  certain  amount  of  computational 
geometry  was  required  near  the  boundary  in  every 
time  step. 

We  consider  the  inviscid  Euler  equations  of  gas 
dynamics  in  two  space  dimensions 

“,+/(«)„ +  g(u)f=  0.  (I) 


where 


‘  P  " 
P  ut 
pu2 
P  E 


f(u) 


pu, 

pu]+p 

pu,u2 

u\(pE  +P) 


g(u)  = 


pu, 
pu,u2 
pu(  +  p 
u,(p£  +p) 


(2) 


Here  (m,,u2)  represents  the  velocity,  E  is  the  total 
energy  per  unit  mass,  and  p  is  the  pressure  which  is 
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related  to  the  other  variables  by  the  equation  of  state. 
We  assume  a  y-law  gas,  so  that 

P  =(Y-  0(p£  -?P(“T  +  Mj))-  (3) 

At  a  solid  wall  boundary  we  require  that  the 
component  of  velocity  normal  to  the  wall  be  zero. 
In  one  space  dimension  the  system  reduces  to 

«,+/(«),=  0.  (4) 

where  «=(p,  pt\  p£)  and  J\u)  =  (pr.  pr:  +  p,  r(p£ 
+  p)),  with  v  =«,  the  velocity.  The  boundary  con¬ 
ditions  become  v  =  0  at  a  solid  wall. 

2.  A  ONE-DIMENSIONAL  EXAMPLE 

In  order  to  illustrate  this  approach  we  begin  with  a 
one-dimensional  mode!  problem,  the  one-dimensional 
Euler  equations  for  x  >  0  with  a  solid  wall  at  x  =  0. 
We  take  a  grid  with  cell  interfaces  at  the  points 

.t0  =  0 

•v  ,=*■ 

xt  =  A '  +  jh  for  j  =  2,1 ... . 

Here  A  is  a  uniform  grid  spacing  and  A'  <  A.  The  grid 
is  uniform  except  for  one  small  cell  near  the  boundary 
(see  Fig.  1).  We  use  a  conservative  method  in  the 
form 

uy‘  =  u;-^[F;,,-F;].  j  =  o.i...  (5) 

h, 

Here  A,  is  the  width  of  the  yth  cell,  so  in  our  case  we 
have  A0  =  A'  and  A;  =  A  for  j  >  0. 

For  simplicity  we  restrict  our  attention  to 
Godunov’s  method  in  the  regular  portion  of  the  grid, 
although  the  ideas  we  propose  can  be  extended  to 
higher  order  methods  as  well.  In  Godunov's  method 
we  take 

F;=/(«*(f/;.,,  u;)),  (6) 

where  u*(uL,u*)  represents  the  solution  to  the 
Riemann  problem  with  left  and  right  states  uL  and 
uR,  evaluated  along  x  t=  0.  Although  a  rigorous 
stability  proof  is  not  available  for  systems  of 
equations,  in  practice  this  method  is  always  stable 
provided  the  CFL  condition 


Fig  I .  One-dimensional  grid  with  one  irregular  cell  adjacent 
to  the  wall. 


is  satisfied,  where  is  the  maximum  wave  speed. 
We  will  assume  that  our  time  step  k  is  chosen  so  that 
the  condition  (7)  is  satisfied  relative  to  the  uniform  A. 
We  will  use  the  flux  Eq.  (6)  for  y  =  2,  3  . . . ,  i.e.  at  all 
interfaces  where  the  cell  on  both  sides  is  regular.  Our 
task  is  to  define  fluxes  F"  for  y  =  0.  1  so  that  we 
maintain  stability  (and  accuracy)  with  this  time  step 
even  if  A’  A. 

First  suppose  h'  =  h.  Then  we  can  use  the 
Godunov  flux  [Eq.  (6)]  also  at  y  =  1.  At  the  wall  we 
use  the  well-known  observation  that  the  solution  to 
the  boundary  value  problem  can  be  obtained  by 
ignoring  the  wall  and  extending  the  computational 
domain  to  the  whole  line  —  x  <  x  <  x,  if  we  take 
data  u0(x)  for  v  <  0  equal  to 

p(x.  0)  =  p(  —x.  0) 

r(.v,  0)  =  —  r(  —  x,  0)  for  v  <  0 

p(x,  0)  =p(-x.  0). 

We  will  denote  this  "reflection"  of  the  data  (in  which 
the  velocity  is  negated)  by  the  operator  .rf,  so  that  : 
shorthand  we  can  write 

u(x,  0)  =  Jt(u( -  x,  0))  for  t  <  0. 

With  this  extended  data,  the  solution  continues  to 
satisfy  u(x,  t)  =  J(u(  —  x.  t))  also  for  />0  and  in 
particular  the  boundary  condition  u(0,  /)  =  ()  is 
automatically  satisfied.  This  suggests  that  we  obtain 
a  flux  at  the  wall  by  solving  a  Riemann  problem  with 
left  and  right  states 

uL  =  J(Un)  uR  =  V a 

in  each  time  step  to  obtain 

Flt=J(u*(J(U0).un)). 

(For  brevity  we  will  leave  off  the  superscript  n  in 
general.)  Note  that  the  density  and  energy  component' 
of  this  flux  will  be  zero  since  the  velocity  component 
u*  is  zero  at  the  wall.  There  will  only  be  a  momentum 
flux  at  the  wall  due  to  the  pressure  there,  as  expected 
physically. 

If  A  '  <  A  we  could  attempt  to  use  this  formula  to 
define  F„  but  we  would  find  that  it  is  unstable  unless 
the  CFL  condition 


is  satisfied.  This  will  place  an  unreasonable  restriction 
on  k  if  A'  4,  A. 

This  instability  is  caused  by  the  fact  that  the 
boundary  flux  Fn  is  based  on  the  data  tj,  alone  It 
the  CFL  condition  (8)  is  satisfied,  then  it  is  only  this 
data  that  affects  the  flux  at  the  wall  over  the  time  siep 
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However,  when  Eq.  (8)  is  violated  the  value  £7,  should 
also  affect  the  flux  at  the  wall,  and  ignoring  this  effect 
leads  to  instability. 

In  a  "large  time  step”  approach  we  increase  the 
stencil  of  the  method,  meaning  we  allow  more  data 
points  to  come  into  the  computation  of  each  flux,  and 
hence  retain  stability.  One  way  to  achieve  this  is  by 
a  wave  propagation  approach.  The  solution  of  the 
Riemann  problem  at  each  cell  interface  consists  of 
three  waves  propagating  away  from  the  interface.  If 
Eq.  (8)  is  satisfied  then  these  waves  remain  in  the  cells 
bordering  the  interface  during  the  entire  time  step 
and  hence  affect  the  solution  only  in  these  cells.  If  Eq 
(8)  is  violated  then  the  waves  may  affect  cells  further 
away.  Implementing  Godunov's  method  in  terms  of 
this  wave  propagation  approach,  allowing  waves  to 
affect  more  than  just  the  adjacent  cell,  gives  a  large 
time  step  generalization  that  remains  stable  for  much 
larger  time  steps.'  In  the  present  context  this  allows 
us  to  reduce  h'  without  reducing  the  time  step  k. 
Waves  from  the  boundary  Riemann  problem  cross 
the  interface  at  v,  and  affect  17,  as  well  as  U0.  Waves 
from  the  interface  at  v,  may  reach  the  boundary. 
These  waves  reflect  off  the  boundary  and  the  reflected 
wave  affects  the  value  and  perhaps  also  Ut  if  the 
reflected  wave  reaches  the  cell  interface  at  v,  during 
the  time  step. 

A  more  detailed  description  of  this  procedure  may 
be  found  in  Ref.  3.  A  natural  extension  to  two  space 
dimensions  gives  one  method  to  deal  with  small  cells 
near  the  boundary,  as  described  in  Refs  1-3.  In  one 
dimension  this  works  very  well  but  in  two  dimensions 
occasional  stability  problems  have  still  been  observed 
due  to  multidimensional  effects. 

2. 1 .  The  new  approach 

Our  new  approach  to  the  small  cell  problem  can 
also  be  illustrated  with  the  one-dimensional  problem 
described  above.  We  again  use  the  method  of  Eq.  (5) 
with  Godunov  fluxes  [Eq.  (6)]  for  y  =  2,  3  ... .  For 
j  =  0  and  7  =  1  we  define  fluxes  in  a  similar  manner 
but  with  a  different  choice  of  states  uL  and  uR  in  the 
Riemann  solver.  Recall  that  in  a  naive  attempt  to  use 
Godunov's  method  regardless  of  the  size  h'  of  the 
small  cell  we  would  take  left  and  right  states 
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(9) 

u\  =  L'„  uR  =  L\  . 

(10) 

To  maintain  stability  when  h‘  is  small,  we  need  to 
allow  data  from  additional  grid  cells  to  affect  the 
left  and  right  states  at  each  of  these  interfaces.  Recall 
that  the  method  is  assumed  to  be  stable  with  our 
choice  of  k  and  h  on  the  regular  portion  of  the  grid. 
This  suggests  that  we  would  define  uj  by  taking  the 
average  value  of  U  over  an  interval  of  length  h  to 
the  left  of  the  interface  v,  and  define  uR  by  taking  the 
average  value  of  U  over  an  interval  of  length  h  to  the 
right  of  xr 


For  example,  at  x0  =  0  (the  wall)  we  set 

=  (h'U0  +  (h  -h ')£/,).  (11) 

If  we  view  the  grid  values  as  defining  a  piecewise 
constant  function  with  values  U,  in  the  yth  cell,  then 
Eq.  ( 1 1 )  is  the  average  value  of  this  function  over  the 
interval  0  ^  .x  $  h.  Note  that  if  h‘  =  h  (the  grid  is 
completely  regular)  then  Eq.  (11)  reduces  to  uR  =  i',t 
as  expected  for  Godunov's  method.  Recall  that  in 


Godunov's  method  we  take  =  J(Ua)  =  Jt{uR ) 
to  impose  the  boundary  condition  r(0,  ;)  =  0.  This 
suggests  that  more  generally  we  take 

u,t  =  ^(ttj),  (12) 

where  uR  is  defined  by  Eq.  (11).  We  then  use  the 
Godunov  flux 

F0=f(u*{uk.u$))  (13) 


as  the  flux  at  the  wall.  Using  Eq.  (12)  guarantees  that 
there  will  be  no  flux  of  mass  or  energy  through  the 
wall  and  hence  that  the  method  is  conservative. 

To  define  the  left  and  nght  states  of  .r,  we  again 
construct  intervals  of  length  h  to  either  side  of  this 
point  and  average  the  piecewise  constant  function 
defined  by  U  over  these  intervals.  To  the  nght  of  .v, 
lies  a  regular  cell  of  length  h,  and  so 

V,.  (14) 

To  the  left  of  .v,  an  interval  of  length  h  extends 
beyond  the  wall  (assuming  h'  <  h).  Beyond  the  wall 
we  assume  that  (J  takes  the  value  given  by  Eq.  (12). 
A  weighted  average  of  this  value  and  UQ  gives  u\ 

u[  = -(h  'U0  +  (h  -  h')uk).  (15) 

h 

The  flux  )\  is  then  defined  by 

F,  =/(m*(u[,  uf )).  (16) 

Again,  if  h '  =  h  this  reduces  to  the  standard  Godunov 
flux. 

This  method  remains  stable  even  when  h'  <?  h.  To 
see  why  this  should  be  so.  consider  the  formula  (5)  for 
7=0  where  h,  =  h It  is  the  division  by  h '  that  may 
cause  stability  problems  unless  the  fluxes  F„  and  Fj 
themselves  agree  to  O(h')  as  /r '  — »0.  The  Godunov 
fluxes  based  on  Eqs  (9)  and  (10)  do  not  have  this 
property.  However,  our  proposed  fluxes  (13)  and  (16) 
do  have  this  property,  since  inspection  of  the  formulas 
(II).  (12).  (14)  and  (15)  shows  that  wf  =  +  0{h') 

and  uR  =  +  O(h')  as  h'-*0.  Since  the  flux  function 

f(u*{uL ,  u "))  is  a  Lipschitz  continuous  function  of  uL 
and  uR,  it  follows  that  Fj  —  F„  =  O(h')  as  /T->0  and 
there  is  at  least  a  chance  that  the  method  remains 
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stable  for  arbitrary  h'4h.  Numerical  experiments 
show  that  this  is  indeed  the  case  (although  it  is  poss¬ 
ible  to  contrive  examples,  such  as  a  strong  rarefaction 
wave  originating  at  this  irregularity,  where  the  results 
are  not  very  accurate). 


3.  BOUNDARY  CONDITIONS  IN  TWO  DIMENSIONS 


Turning  now  to  the  two-dimensional  problem,  we 
will  give  a  brief  description  of  how  the  idea  described 
above  extends  to  handle  the  small  cell  problem. 

Consider  the  portion  of  the  boundary  shown  in 
Fig.  2a  and  a  typical  boundary  cell  (i.j).  The  formula 
for  updating  the  value  U ,,  is  the  two-dimensional 
analog  of  Eq.  (5) 

Ul*  ‘  =  V„  —  -j-  [F, * ,,  -  F„  +  G,j  . ,  -  G„  +  #/„]. 

(17) 


The  fluxes  F,  G.  and  H  represent  flux  per  unit  time 
through  the  corresponding  side  of  the  grid  cell  (see 
Fig.  2b)  and  A,t  is  the  area  of  the  cell.  If  any  of  the 
sides  are  missing,  the  corresponding  flux  is  zero. 

On  regular  grid  cells,  //,(  =  0  and  the  fluxes  F  and 
G  might  be  defined  by  an  extension  of  the  Godunov 
method,  setting 


F,,  —  hf(u*(U, .  | 

G,,  =  hg(u*(Ult _ 

Here  u*  represents  the  solution  to  the  appropriate 
Riemann  problem  in  the  v  or  y  direction.  Note  that 
the  fluxes  include  the  factor  h,  the  length  of  each  side, 
to  give  a  flux  per  unit  time  across  the  side. 

It  is  the  denominator  A„  in  Eq.  (17)  that  causes 
trouble  when  the  cell  is  very  small.  We  again  assume 
the  method  is  stable  on  the  regular  portion  of  the 
grid,  where  A(J  =  h2.  To  maintain  stability  we  need  to 
insure  that  our  formulas  for  the  fluxes  cause  the  total 
flux  [the  sum  in  brackets  in  Eq.  (17)]  to  cancel  to 
0(At/)  as  Av-* 0.  This  is  only  possible  if  the  fluxes  are 
computed  via  formulas  that  involve  more  than  just 
the  two  cells  bordering  the  cell  side.  We  take  an 
approach  analogous  to  what  we  described  above  in 
one  dimension. 


Fig.  3.  The  inbox  and  outbox  constructed  from  the  bound¬ 
ary  segment  of  cell  (i.j).  and  the  inbox  for  two  neighboring 
cells. 


sions  the  solid  wall  boundary  condition  requires  that 
the  normal  velocity  at  the  wall  be  equal  to  zero.  If  we 
have  some  value  u\°  representing  the  value  oi  <■  just 
inside  the  wall,  then  we  can  obtain  the  flux  H,,  by 
solving  a  one-dimensional  Riemann  problem  in  the 
direction  normal  to  the  wall,  with  left  and  right  states 

uf;  =  Jf(u^)  u*  =  u™. 

The  reflection  operator  A  is  now  defined  by  negating 
the  normal  velocity  component  while  leaving  the 
tangential  velocity  component  along  with  the  density 
and  pressure  unchanged.  The  resulting  Godunov  flux 
is  used  for  H,r 

We  obtain  by  a  procedure  analogous  to  that 
of  the  one-dimensional  example.  We  construct  a  box 
extending  a  distance  h  away  from  the  wall  as  shown 
in  Fig.  3.  The  box  extending  into  the  computational 
domain  is  called  inbox(/,y).  The  mirror  image  box 
outside  the  domain  is  called  outbox(/'.y).  We  obtain 
the  value  u !°  by  viewing  the  given  data  U  as  defining 
a  piecewise  constant  function,  constant  in  each  grid 
cell,  and  setting  to  be  the  average  value  of  this 
function  over  the  region  inbox(i.y).  In  Fig.  3  inbo x(i.y') 
would  contain  an  area-weighted  average  of  two  grid 
values  while  the  value  for  inbox(i  —  I  .j)  is  based  on 
four  grid  values.  We  think  of  the  outbox  as  containing 
the  value  m“ui  =  Jt(u'"). 

To  find  the  weights  needed  to  compute  u'°  we  must 
compute  the  intersection  of  the  inbox  with  each 
nearby  cell.  This  is  easily  accomplished  with  standard 
computational  geometry  routines.  Note  that  for  a 
given  geometry  and  grid  these  weights  need  only  be 
computed  once  at  the  beginning  of  the  computation 
They  need  not  be  recomputed  in  each  time  step. 


3.1.  Boundary  fluxes 

We  begin  by  considering  the  boundary  segment, 
where  we  must  compute  the  flux  H„.  In  two  dimen- 


(b) 


Fi*'i 


Fig.  2.  (a)  The  Cartesian  grid  near  the  boundary  lb)  Blow¬ 
up  of  cell  (i.j)  showing  the  location  of  fluxes 


3.2.  Fluxes  at  other  sides 

We  now  consider  the  fluxes  F  and  G  along  other 
sides  of  this  cell.  These  are  all  computed  by  similar 
procedures,  so  to  be  specific  we  will  consider  the 
computation  of  F„.  the  flux  on  the  left  side  of  this  cell. 

To  compute  Fu  we  solve  two  Riemann  problems, 
one  in  some  direction  %  with  data  uf .  uf  and  the  other 
in  the  orthogonal  direction  q  with  data  u2.  uj. 
The  choice  of  these  directions  and  the  data  will  be 
discussed  in  a  moment.  First  we  explain  how  these 
Riemann  problem  solutions  are  computed  and  used 
to  define  F„. 
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Fig.  4.  A  vertical  cell  interface  and  the  and  n-directions 

Figure  4  shows  a  typical  vertical  cell  interface 
and  two  orthogonal  directions  \  and  13.  Let  9  be  the 
angle  that  E,  is  rotated  from  the  x -direction  (9  <  0  in 
this  example).  Suppose  we  solve  a  one-dimensional 
Riemann  problem  in  the  ^-direction  with  left  and 
right  states  u£,  ut  to  obtain  the  flux  per  unit  length 
per  unit  time  in  the  ^-direction.  (To  do  this  we  rotate 
the  velocity  components  of  ub.  into  t)  velocity 
components,  solve  the  one-dimensional  Riemann 
problem,  and  then  rotate  the  resulting  flux  /  back 
to  x-y  velocity  components.)  Call  his  resulting 
flux/,. 

Similarly,  solving  a  one-dimensional  Riemann 
problem  in  the  rj-direction  with  left  and  right  slates 
u%,  u*  gives  /n,  the  flux  per  unit  length  per  unit  time 
in  the  r|-direction.  The  total  flux  across  the  vertical 
segment  of  length  h '  is  then 

F  -  h  '(/=  cos  0  -  j\  sin  9).  (19) 

This  is  the  value  we  use  for  the  flux  F,r 

This  same  approach  has  been  used  by  others 
(see  Refs  8-10)  to  define  multidimensional  upwind 
methods.  In  these  methods  the  directions  %  and  13 
are  chosen  based  on  the  local  flow  in  an  attempt  to 
use  physically  meaningful  directions  in  place  of  the 
artificial  coordinate  directions.  For  example,  the 
direction  of  the  velocity  or  the  pressure  gradient 
might  be  used  to  define  t,.  In  our  application  we  are 
only  considering  cells  adjacent  to  the  boundary  and 
the  relevant  directions  are  the  directions  tangential 
and  normal  to  the  wall.  We  choose  ^  to  be  the 
direction  tangential  to  the  wall  in  one  of  the  two  cells 
bordering  this  interface.  Since  our  primary  concern 
is  to  maintain  stability  in  very  small  cells,  we  choose 
the  smaller  of  the  two  adjacent  cells  to  define  this 
direction.  This  will  lead  to  cancellation  of  fluxes  in 
tiny  cells  in  the  same  manner  as  previously  seen  in  the 
one-dimensional  example.  The  r|-direction  is  normal 
to  the  ^-direction. 

3.3.  Tangential  boxes 

We  must  still  specify  the  data  for  these  tangential 
and  normal  Riemann  problems.  We  first  consider  the 
tangential  problem.  We  use  an  approach  similar  to 
the  specification  of  data  in  an  inbox  described  above. 
From  the  interface  we  construct  boxes  that  extend  a 
distance  h  in  the  ^-direction.  Figure  5a  shows  an 
example.  The  data  u f.  u f  are  obtained  by  an  area- 
weighted  average  of  the  values  in  each  cell  the  bc  \ 


Fig.  5.  (a)  Tangential  boxes  constructed  from  the  cell 
interface,  (b)  Normal  boxes  fiu„.  he  cell  interface  and 
outboxl  1  -  t.y) 

overlaps.  In  our  current  implementation  we  assume 
the  wall  is  convex,  so  that  these  boxes  lie  entirely 
wtihin  the  computational  domain.  Each  box  overlaps 
at  most  two  grid  cells  and  the  weights  are  easily 
calculated.  Since  the  directions  -  and  13  and  the 
resulting  boxes  depend  only  on  the  geometry,  not  on 
the  flow  variables,  these  weights  can  again  be  calcu¬ 
lated  once  and  for  all  as  a  preprocessing  step. 

3.4.  Normal  boxes 

Figure  5b  shows  the  normal  boxes  in  the  13-direc¬ 
tion.  The  box  in  the  outward  direction  does  not  hit 
the  boundary  and  overlaps  at  most  two  regula  v  ‘s. 
so  u *  is  calculated  as  an  area-weighted  averag.  of 
these  cell  values.  The  other  box  may  extend  beyond 
the  boundary.  If  so,  the  portion  lying  outside  the 
computational  domain  lies  in  one  or  more  outboxes. 
the  artificial  cells  created  in  the  process  of  computing 
the  boundary  flux  H,t  described  above.  Figure  5b 
shows  a  simple  example  where  the  normal  box 
intersects  only  one  cell  (/  -  1 .7)  and  outbox(i  -  l.y). 
More  generally  the  normal  box  might  intersect  two 
cehs  and  their  outboxes.  as  happens  for  example 
when  we  compute  the  flux  F, , , ,  which  involves  cells 
(t\j)  and  (/.  /  -  1 ).  Moreover  the  two  outboxes  will  in 
general  overlap  due  to  the  convexity  of  the  region. 
We  again  use  area-weighted  averaging  over  the  Tour 
cells  in  question,  weighting  the  values  U  ,  C, ,  , .  . 

w”“‘  |  by  the  areas  of  intersection  and  then  dividing  by 
the  sum  of  all  these  areas. 

3.0.  Cancellation 

Although  we  will  rot  present  the  details  here,  it  can 
be  shown  that  this  wav  of  defining  fluxes  leads  to  the 
desired  cancellation  ol  fluxes  in  very  small  cells.  1 .  » 
values  of  ub  computed  at  each  of  the  three  sides  of  a 
very  small  triangular  cell  are  nearly  the  same  because 
of  inr  construction.  Th"y  differ  by  only  0(.4, ,)  as 
A„-  0  The  same  is  true  of  each  of  the  other  values 
u ?,  u[ ;.  »*  and  so  by  Lipschitz  continuity  of  the  fluxes 
F,  G  and  H  we  obtain  the  required  cancellation. 
Numerical  results  show  stability  even  when  4„  is 
many  orders  of  magnitude  less  than  /r 

3.6.  Higher  order  methods 

The  method  (17)  with  fluxes  ( 18)  is  only  first  order 
accurate  and  ii  highly  dissipative.  In  our  previous 
work  we  used  the  wave  propagation  boundary  con¬ 
ditions  together  with  a  high  resolution  method  away 
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from  the  boundary  and  obtained  reasonable  results 
(see  Ref.  1).  The  new  boundary  conditions  can  also 
be  applied  in  conjunction  with  a  high  resolution 
method  and  give  similar  results.  Moreover,  with  our 
new  formulation  it  appears  to  be  easier  to  improve 
the  accuracy  of  the  boundary  conditions,  allowing 
us  to  obtain  higher  order  accuracy  overall.  The  main 
idea  is  to  introduce  slopes  in  each  cell  and  use 
piecewise  linear  approximations  in  place  of  p.ccewise 
constants  to  define  the  fluxes.  Near  the  boundary  we 
can  easily  estimate  slopes  in  the  tangential  direction 
along  the  wall  by  differencing  values  in  the  inboxes 
that  we  have  defined  above. 

These  improvements  are  still  being  investigated 
and  will  be  reported  in  detail  elsewhere.  Here  we  will 
only  compare  results  obtained  with  the  method  as  we 
have  described  it  and  results  obtained  using  the  same 
interior  method  with  the  wave  propagation  boundary 
conditions  described  in  earlier  papers. 

4.  NUMERICAL  results 

We  show  one  representative  test  case,  a  supersonic 
shock  going  around  an  expansion  corner.  We  also 
show  the  steady  state  solution  obtained  at  la^ge  times. 
The  exact  rarefaction  wave  solution  is  a  simple  wave 
and  can  be  computed  following  Sec.  6.17  of 
Whitham,"  for  example. 

The  geometry  we  use  is  tne  rectangle  [0,  1.32]  x 
[0.  0.8]  with  a  solid  wall  at 

f0.3  v^O.l 

y  =  0.3(1 -(.r- 0.1 );)  0. 1  <  .v  <  0.7 

[o.  192  —  0.36(.v  -0.7)  0.7  sg  x  $  1.32. 


The  initial  conditions  consist  of  a  Mach  2.31  shock 
at  x  =  0.06  with  left  and  right  states 

pL  =  5.1432,  u[  =  2.045 1 1.  ui  =  0. 
pL  =  9. 04545 
and 

p*=l.4.  wf  =  0.  m?  =  0.  pR  =  1.0. 

We  take  h  =  0.02  (66  x  40  grid)  and  a  time  step 
k  =  0.002.  This  corresponds  to  a  Courant  number  of 
roughly  0.37  relative  to  the  regular  cells  with  area  > 
For  the  crude  form  of  Godunov's  method  used  here, 
the  stability  restriction  requires  a  Courant  number  of 
less  than  0.5.  The  smallest  cells  near  the  boundary 
have  an  area  of  roughly  10” 3  h :. 

Figure  6  shows  numerical  results  at  time  t  =  0.4.  as 
the  shock  is  rounding  the  corner.  Results  obtained 
with  the  wave  propagation  boundary  conditions  are 
shown  in  Fig.  6a,  while  Fig.  6b  shows  the  results 
obtained  with  our  new  approach.  These  results  are 
very  similar.  Slight  discrepancies  can  be  seen  near  the 
wall  just  around  the  shock.  For  this  problem,  both 
sets  of  boundary  conditions  worked  well.  We  have 
also  performed  tests  on  other  problems  where  the 
wave  propagation  method  shows  instabilities  and 
have  observed  no  such  difficulties  with  the  new 
method. 

Figure  7  shows  the  steady  state  results  obtained 
after  many  iterations  of  the  time  dependent  code  (no 
attempt  has  been  made  so  far  to  accelerate  converg¬ 
ence  for  steady  state  solutions).  We  only  show  the 
results  with  our  new  boundary  conditions.  The  wave 
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Fig.  6  Shock  propagation  results  at  r=0.4.  (a)  Using  the  wave  propagation  boundary  conditions. 

(b)  Using  the  new  boundary  conditions 


Fig  1.  Steady  state  results,  (a)  Pressure  contours,  (b)  Pressure  along  the  wall  The  solid  line  is  the  exa'"t 
olution.  +  indicate'  the  numerical  solution. 
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propagation  boundary  conditions  give  very  similar 
results.  We  use  a  coarser  grid  than  in  the  previous 
example  (h  =  0.04)  in  order  to  demonstrate  that  we 
achieve  reasonable  accuracy  along  the  boundary  even 
with  a  relatively  coarse  piecewise  linear  representation 
of  the  boundary.  We  also  use  a  larger  computational 
domain,  [0.  2]  x  [0,  1.6],  to  minimize  the  impact  of  the 
far-field  boundaries.  The  true  solution  is  a  rarefaction 
wave  originating  from  the  portion  of  the  boundary 
with  nonzero  curvature.  In  the  exact  solution  the 
contour  lines  would  be  straight  lines.  Our  results  are 
contaminated  by  effects  from  the  far-field  boundary. 

Near  the  solid  wall  the  contour  lines  appear  to 
show  a  boundary  layer.  This  is  an  artifact  of  the 
graphics  routine,  which  assumes  the  data  are  on 
a  uniform  grid  at  cell  centers.  Our  data  near  the 
boundary  should  be  viewed  as  an  approximation  to 
the  pointwise  value  at  the  center  of  mass  of  the 
irregular  cell,  not  at  the  center  of  the  full  Cartesian 
cell. 

In  order  to  examine  the  accuracy  at  the  wall. 
Fig.  7b  shows  plots  of  the  pressure  along  the  bound¬ 
ary.  plotted  against  arclength.  To  obtain  a  boundary 
pressure,  the  cell  value  L'„  and  the  reflected  value 
3?(  are  used  to  solve  the  one-dimensionel 
Riemann  problem  normal  to  the  wall  in  each  irregular 
cell.  The  resulting  pressure  p*  is  used  as  the  boundary 
pressure.  Figure  7b  shows  these  results  and  also  the 
exact  solution  (to  machine  precision)  calculated  using 
the  theory  of  Ref.  1 1. 

In  more  complicated  computations  we  use  adaptive 
grid  refinement  to  obtain  high  resolution  results  with 
minimal  effort.  The  boundary  conditions  described 
here  can  also  be  used  in  conjunction  with  the  adaptive 
Cartesian  grid  code  described  in  Refs  1  and  2. 
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The  aim  of  this  work  is  the  development  of  an  automatic,  adaptive  mesh  refinement  strategy 
for  solving  hyperbolic  conservation  laws  in  two  dimensions.  There  are  two  main  difficulties  in 
doing  this  The  first  problem  is  due  to  the  presence  of  discontinuities  in  the  solution  and  the 
effect  on  them  of  discontinuities  in  the  mesh.  The  second  problem  is  how  to  organize  the  algo¬ 
rithm  to  minimize  memory  and  CPU  overhead.  This  is  an  important  consideration  and  will 
continue  to  be  important  as  more  sophisticated  algorithms  that  use  data  structures  other  than 
arrays  are  developed  for  use  on  vector  and  parallel  computers.  1  iw  Academic  Press.  Inc 


1.  Introduction 

In  this  paper,  we  present  computations  that  use  adaptive  mesh  refinement 
to  solve  multidimensional,  time  dependent  shock  hydrodynamics  problems. 
Complicated  structures  such  as  multiple  Mach  reflections  arise  in  these  problems. 
Adaptive  techniques  are  essential  for  our  computations  in  order  to  adequately 
resolve  features  in  the  solution  within  today's  computer  limitations. 

Our  starting  point  will  be  the  algorithms  in  [6]  for  adaptive  mesh  refinement  for 
hyperbolic  equations  on  rectangular  grids.  In  this  approach,  the  refined  regions 
consist  of  a  small  number  of  rectangular  grid  patches  with  finer  mesh  spacing  than 
the  underlying  global  coarse  grid.  These  rectangular  subgrids  contain  points  where 
the  error  in  the  coarser  grid  solution  is  too  high,  and  other  points  as  well.  We  use 
rectangular  subgrids  so  that  we  can  use  integration  methods  for  rectangular  grids 
whose  convergence  properties  are  well  understood.  These  methods  can  be  made 
quite  efficient  on  vector  and  parallel  computers.  In  addition,  rectangular  grids  have 
a  simple  user  interface.  We  can  use  the  same  integrator  on  fine  and  coarse  grids.  By 
separating  the  integrator  from  the  adaptive  strategy,  an  off-the-shelf  integrator  can 
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be  used  without  modification.  This  eliminates  much  of  the  problem  specific  work  in 
doing  adaptive  calculations. 

The  present  work  differs  from  that  in  [6]  in  several  respects.  The  main  one  is 
that  we  are  computing  unsteady  flows  with  shocks,  so  that  maintaining  global  con¬ 
servation  form  is  a  primary  consideration.  The  second  difference  is  that  the  nested 
refinements  we  use  have  boundaries  coinciding  with  the  grid  lines  of  the  underhing 
coarse  mesh.  This  greatly  simplifies  the  maintenance  of  conservation  over  the 
former  approach,  where  the  refined  subgrids  were  allowed  to  be  rotated  with 
respect  to  the  coarse  grid.  Third,  great  care  was  taken  to  obtain  an  efficient 
implementation  on  a  supercomputer.  The  main  difficulty  with  adaptive  method-  i- 
the  need  for  data  structures  not  usually  found  in  numerical  software.  We  felt  the 
program  complexity  was  high  enough  to  justify  the  effort  of  devising  as  general  and 
automatic  an  approach  as  possible. 

Earlier  work  along  the  lines  of  the  present  work  was  done  by  ["]  in  one  dimen¬ 
sion.  and  [14]  for  scalar  problems  in  two  dimensions.  [19]  have  also  computed 
transonic  flow  in  two  dimensions  with  grid  embedding.  However,  in  the  latter  ;w. 
approaches,  the  grids  were  not  restricted  to  rectangles.  The  data  structure-,  and 
therefore  the  efficiency  of  such  an  approach,  are  quite  different.  Our  method  of 
adaptivity  through  grid  refinement  is  in  contrast  to  methods  that  adapt  the  grid  K\ 
moving  grid  lines  into  one  region,  leaving  a  coarser  region  somew  here  else  [ !  v  • 
15.  13.  20.  8].  Such  methods  try  to  get  the  most  accurate  solution  for  a  lived  _ 
whereas  our  approach  tries  to  attain  a  fixed  accuracy  for  a  minimum  co-t  B. 
approaches  have  their  advantages  and  disadvantages.  The  so-called  moving  g-d 
point  methods  often  have  trouble  maintaining  a  smooth  grid.  Regularity  term-  m. 
penalty  functions  used  to  regularize  the  grid  add  overhead  and  reduce  the  -imp:  „ .: . 
of  these  methods.  Local  grid  refinement,  on  the  other  hand,  has  the  drawb.uk  • 
needing  special  equations  at  grid  interfaces.  In  a  method  where  a  fixed  nunr.-: 
grid  points  are  used  during  a  computation,  the  user  must  initially  guess  at  wh.r  w  ; 
be  an  adequate  number  of  points  to  resolve  features  in  the  solution  that  may 
later  With  local  grid  refinement,  grid  points  are  added  or  removed  as  ncco— urv 

In  the  numerical  experiments  shown  below,  we  have  combined  thi-  adam 
mesh  refinement  strategy  with  the  high  resolution  difference  scheme  of  [I"  • 

develop  an  almost  automatic  software  tool  for  solving  gas  dynamics  proniem- 
two  space  dimensions.  A  reasonable  question  is.  why  is  an  adaptive  method  ne.  tv.: 
given  that  the  difference  scheme  used,  a  second-order  Godunov -type  me' - 
already  has  quite  high  resolution'1  Conversely,  if  adaptive  methods  are  u-ed  .-. 
such  complicated  and  expensive  difference  schemes  really  necessary ?  The  an-sc  - 
that  both  components  are  necessary  to  obtain  well-resolved  results  tor  -  ,u- 
hydrodynamics.  It  has  been  demonstrated  [22]  that  the  more  complu..:--: 
Godunov-type  schemes  give  more  resolution  per  computational  dollar  than  -imp  v 
schemes  -uch  as  Lax-Wendroff.  Given  that  a  high  quality  scheme  is  neev--.;:  . 
adaptive  mesh  refinement  can  then  concentrate  the  computational  effort  in  regi<m- 
where  it  is  most  useful.  Since  Godunov-type  methods  are  mere  expensive  tli.m 
simple  schemes,  the  computational  savings  of  selective  refinement  can  be  -ub-tan- 
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tial.  Some  of  the  computations  presented  here  could  not  have  been  done  reasonably 
without  the  use  of  an  adaptive  solver. 

In  the  sections  that  follow,  we  describe  the  adaptive  mesh  refinement  (AMR) 
algorithm  for  integrating  a  general  hyperbolic  system  of  conservation  laws 

u,  +  /(u)t  +g{u),  =0  on  Dc.R 2 
Bu  =  b  on  cD. 

Our  numerical  examples  involve  the  Euler  equations  for  gas  dynamics,  where 

(pu 
pu-+p 
puv 

puE  +  up 

and 

P  =  (y -i)(p£- p(“  jj. 

Although  our  work  to  date  is  in  two  space  dimensions,  all  the  algorithms  extend 
to  three  dimensions,  and  in  fact  it  seems  possible  to  implement  a  general  code 
where  the  number  of  dimensions  is  input. 

In  the  next  sections  we  describe  in  detail  the  adaptive  mesh  refinement  algorithm 
and  its  implementation.  We  give  enough  detail  for  users  interested  in  modifying  the 
algorithm,  using  our  code,  or  writing  their  own.  First,  we  discuss  the  structures  that 
define  our  grid  hierarchy.  Next,  we  describe  the  integration  scheme  for  such  a 
(static)  grid  hierarchy.  Third,  the  grid  generation  and  error  estimation  procedures 
used  to  generate  the  grid  hierarchy  itself  are  presented.  Our  error  estimation  proce¬ 
dure  is  theoretically  justifiable  only  for  smooth  solutions.  We  discuss  variations  of 
it  that  may  prove  useful  for  problems  with  shocks.  In  the  last  section  we  present 
numerical  experiments  along  with  a  detailed  timing  analysis  of  the  runs.  This 
program  is  being  used  to  study  Mach  reflection  in  two  dimensions  with  resolution 
not  previously  possible.  New  results,  a  triple  Mach  stem  configuration  at  low 
have  been  observed. 


g(M)  = 


2.  Grid  Description 

AMR  is  based  on  using  a  sequence  of  nested,  logically  rectangular  meshes  on 
which  the  pde  is  discretized.  In  this  work,  we  require  the  domain  D  to  be  a  finite 
union  of  rectangles  whose  sides  lie  in  the  coordinate  directions.  We  assume  here 
that  all  the  meshes  are  physically  rectangular  as  well,  although  this  is  not  essential 
The  method  discussed  here  can  be  implemented  on  a  general  quadrilateral  mesh. 
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(See.  for  example,  [5]).  We  define  a  sequence  of  levels  /  =  1 . /m„.  A  grid  G, .  has: 

mesh  spacing  A,,  with  level  !  coarsest,  and  define 


With  an  abuse  of  terminology  describing  a  grid  and  the  domain  it  covers,  ne  have 
G,  =  G  ,  =  D.  the  problem  domain.  If  there  are  several  grids  at  level  I.  the  grid 
lines  must  "align"  with  each  other;  that  is.  each  grid  is  a  subset  of  a  rectangular 
discretization  of  the  whole  space. 

We  may  often  have  overlapping  grids  at  the  same  level,  so  that  G  -  G  .  *0. 
/  ft  k  However,  we  require  that  the  discrete  solution  be  independent  of  how  G  is 
decomposed  into  rectangles. 

Grids  at  different  levels  in  the  grid  hierarchy  must  be  "properly  nested."  This 
means 

(i  i  a  fine  grid  starts  and  ends  at  the  corner  of  a  cell  in  the  next  coarser  grid, 
nil  There  must  be  at  least  one  le'e!  /—  I  cell  in  some  level  /-  I  grid  separat¬ 
ing  a  grid  cell  at  level  /  from  a  cell  at  level  t  -  2.  in  the  north,  south,  east,  and  west 
directions,  unless  the  cell  abuts  the  physical  boundary  of  the  domain. 

Note  that  this  proper  nesting  is  not  as  stringent  as  requiring  a  fine  grid  to  he 
contained  in  only  one  coarser  level  grid.  For  example,  in  Fig.  2.1.  there  is  one  grid 
at  level  3.  G,  ,.  Every  grid  point  in  G,  .  is  contained  in  one  of  the  two  grids  jt 
level  2.  G;  .  or  G: ;. 

Grids  will  be  refined  in  time  as  well  as  space,  by  the  same  mesh  refinement  ratio 
r.  where  r  =  J.v  ,  J.v.  Thus. 


Jt,  _  Jt  _  _  J tj_ 

J.v  .J.v  J.v, 

and  so  the  same  explicit  difference  scheme  is  stable  on  all  grids.  This  means  more 


Ik.  J  I  Grid  C,  .  spans  two  coarser  grids  but  is  properlv  level  nested 
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time  steps  are  taken  on  the  finer  grids  than  on  the  coarser  grids.  This  is  a 
reasonable  requirement  from  the  point  of  view  of  accuracy,  since  for  many  dif¬ 
ference  schemes,  the  leading  terms  in  the  spatial  and  temporal  truncation  error  are 
of  the  same  order,  fn  addition,  the  smaller  time  step  of  the  fine  grid  is  not  imposed 
globally.  In  this  implementation  we  only  allow  an  even  refinement  ratio.  This 
simplifies  the  error  estimation  procedure  described  later,  by  avoiding  the  need  of 
distinguish  between  an  odd  and  even  number  of  grid  points  in  a  grid. 

At  discrete  times  the  grid  hierarchy  may  be  modified.  The  finest  grids  need  to  be 
changed  (“moved,"  deleted  if  necessary)  most  often.  When  grids  at  level  /  are 
changed,  all  finer  level  grids  are  changed  as  well,  but  the  coarser  grids  may  remain 
fixed.  New  grids  at  level  l  may  replace  the  old  ones,  but  they  are  still  subject  to  the 
same  "proper  nesting"  requirement. 

A  point  (  v.  i  )6  D  may  be  contained  in  several  grids.  The  solution  u[x,  y)  will  be 
taken  from  the  finest  grid  containing  the  point.  If  there  are  several  equally  fine  grids 
containing  the  point,  any  fine  grid  value  will  suffice,  since  the  solution  on  the 
intersection  of  overlapping  fine  grids  will  be  identical. 


3.  Integration  Algorithm 

AMR  assumes  there  is  a  basic,  underlying,  conservative,  explicit  finite  difference 
scheme  of  the  form 


~  -  1  Z.J  -  F,  - 1  :.y>  -  1  :  “  i ;). 


(1 


The  values  u,  ,  are  cell-centered  quantities.  Each  cell  is  defined  by  its  four  corner 
grid  points.  If  there  are  no  refined  regions,  then  Eq.  ( 1 ),  augmented  by  the  dis¬ 
cretized  physical  boundary  conditions,  defines  the  time  evolution  on  a  single  grid. 

With  multiple  grids,  each  grid  is  separately  defined  and  has  its  own  solution 
vector,  so  that  a  grid  can  be  advanced  independently  of  other  grids,  except  for  the 
determination  of  its  boundary  values  (see  Section  4).  The  integration  steps  on  dif¬ 
ferent  grids  are  interleaved,  so  that  before  advancing  a  grid  to  time  r  +  Jr,  all  the 
finer  level  grids  have  been  integrated  to  time  r.  Scheme  (1 )  is  still  initially  applied 
on  every  grid  at  every  level,  but  the  results  will  need  to  be  modified  in  case 

(i)  the  cell  is  overiayed  by  a  finer  level  grid:  or 

(ii )  the  cell  abuts  a  fine  grid  interface  but  is  not  itself  covered  by  any  fine  grid. 

In  case  til.  the  coarse  grid  value  at  level  /  -  1  is  defined  to  be  the  conservative 
average  of  the  fine  grid  values  at  level  /  that  comprise  the  coarse  cell.  After  every 
coarse  integration  step,  the  coarse  grid  value  is  simply  replaced  by  this  conservative 
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average,  and  the  value  originally  calculated  using  ( 1 )  is  discarded.  For  a  refinement 
ratio  of  r,  we  define 


u 


coarse 

/•/ 


r  -  l 

I 


n 

k  ~p 


V 


where  the  indices  refer  to  the  example  in  Fig.  3.1.  We  could  define  a  coarse  cell  u 
at  multiples  of  the  fine  time  step  in  the  same  way.  but  this  is  not  necessary  This 
is  equivalent  (within  roundoff  error)  to  redefining  the  coarse  fluxes  around  the 
overlayed  coarse  grid  point  to  be  the  sum  over  the  fine  time  steps  of  all  fine  grid 
fluxes  calculated  on  any  boundary  segment  for  that  cell.  However,  this  implementa¬ 
tion  would  use  extra  storage  to  save  the  fine  grid  fluxes.  By  updating  the  solution 
values  themselves,  no  extra  flux  storage  is  needed. 

In  case  (iil.  the  difference  scheme  ( 1 )  itself  that  is  applied  to  the  coarse  cell  must 
be  modified.  According  to  1 1 ).  the  fine  grid  abutting  the  coarse  cell  has  no  effect. 
However,  for  the  difference  scheme  to  be  conservative  on  this  grid  hierarchs,  the 
fluxes  into  the  fine  grid  across  a  coarse,  fine  cell  boundary  must  equal  the  flux  out 
of  the  coarse  cell.  (This  conservative  procedure  has  been  discussed  by  [["].  A  fuller 
discussion  of  conservation  at  grid  interfaces  is  in  [4].)  We  use  this  to  redefine  the 
coarse  grid  flux  in  case  (ii).  For  example,  in  Fig.  3.2,  the  difference  scheme  at  point 
i,  j  should  be 

*6.  i  H  +  A  f  coarse  1 

1  '  -  1  '  -  1 

F,.;  V  V  Fk „  .  .rlf  -  q  J/.-in4i 


=  U,  ,(t)~ 


Atc 


J.x 


A  f  coarse 

Jv 


[G,.  ,  -  I  ;H)  -  G,  ,  ;(f)]. 


1 2  i 


where  J.v  and  J  v  are  coarse  spatial  step  sizes.  The  double  sum  is  due  to  the  refine¬ 
ment  in  time:  for  a  refinement  ratio  r.  there  are  r  times  as  many  steps  taken  on  the 
fine  grid  as  the  coarse  grid.  If  the  cell  to  the  north  of  (i,j)  were  also  refined,  the 
flux  G.  ,  . ,  :  would  be  replaced  by  the  sum  of  fine  fluxes  as  well. 


Fig.  3  1  The  coarse  cell  value  is  replaced  by  the  average  of  all  the  fine  grid  points  in  that  cell 


70 


BERGER  AND  COLELLA 


LOCAL  ADAPTIVE  MESH  REFINEMENT 


71 


The  boundary  fluxes  SF  are  stored  in  a  vector  associated  with  every  fine  grid.  In 
the  initialization  step  (3),  there  may  be  several  coarse  grids  that  set  SF.  Since  all 
fluxes  calculated  at  a  given  edge  and  level  are  identical  (up  to  roundoff  error!  and 
are  independent  of  the  particular  grid  on  which  they  are  calculated,  we  simply  use 
the  last  value  assigned.  At  the  end  of  a  time  step,  we  may  have  several  fine  grids 
available  to  update  a  given  coarse  cell  edge,  since  overlapping  grids  are  permitted. 
For  this  reason,  we  use  a  matrix  to  indicate  the  edges  of  a  coarse  cell  that  have 
already  been  updated  and  only  perform  the  update  once  for  each  edge.  As  before, 
it  does  not  matter  which  fine  grid  actually  performs  the  update  for  any  given  edge, 
so  the  result  is  independent  of  the  order  in  which  the  fine  grid  list  is  traversed.  Thi> 
modification  is  a  negligible  amount  of  work,  taking  approximately  0.3%  of  a 
typical  run  time.  On  machines  with  a  scatter  gather  operation,  this  should  proceed 
even  faster. 

We  emphasize  that  this  work  is  done  as  a  "fix  up“  step  after  each  grid  is  updated 
using  scheme  ( 1 ).  In  this  way.  the  integrator  can  be  separated  from  the  additional 
work  which  is  needed  because  of  the  grid  hierarchy.  A  new  difference  scheme  can 
be  substituted  by  a  user  unfamiliar  with  and  not  interested  in  the  inner  working' 
of  the  AMR  program 


4.  Boundary  Conditions 

A  discussion  of  boundary  conditions  completes  the  description  of  the  integration 
procedure  on  a  multiple  grid  hierarchy.  Let  the  interior  integration  scheme  have  ; 
stencil  w  hich  is  centered  in  space,  with  d  points  to  each  side.  To  compute  the  new 
time  step.  AMR  provides  solution  values  at  the  old  time  step  on  a  border  oi  cell' 
of  width  d  intersected  with  the  physical  domain.  The  user  must  supply  the  code  : 
compute  any  additional  information  needed  to  implement  the  boundary  condim-m 
For  example,  if  boundary  conditions  are  imposed  by  extrapolation,  the  user  w.>Jd 
provide  the  extrapolated  values  for  points  outside  the  domain. 

For  a  grid  at  level  /.  the  bordering  cell  values  are  provided  using  value'  :r 
adjacent  level  /  grids  where  they  are  available;  otherwise,  the  AMR  algorithm  ->m 
putes  boundary  values  using  bilinear  interpolation  from  coarser  level  >oi 
values.  If  necessary,  we  also  interpolate  linearly  in  time. 

It  may  happen  that  a  point  i.v.  y)  is  inside  the  domain  D.  but  one  or  more 
rounding  coarse  grid  points  needed  for  the  bilinear  interpolation  are  outsid:  \ 
before,  we  assume  there  is  a  user-supplied  routine  that  can  provide  exterior  v.  mw 
grid  points  given  some  interior  points. 

Our  implementation  partitions  the  required  border  cells  at  level  /  into  rectangula: 
boundary  patches.  For  each  rectangular  piece  we: 

in  find  solution  values  from  level  /-  1  grids  on  a  slightly  larger  rectangular 
piece  enclo>ing  the  border  cells; 

in  i  linearly  interpolate  for  the  border  values: 
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(iii)  if  there  are  fine  grids  at  level  /  that  could  supply  some  values  (say  an 
adjacent  fine  grid),  overwrite  the  linearly  interpolated  values  from  step  (ii). 

In  step  (i),  most  of  these  coarser  level  values  are  found  by  intersecting  the  rec¬ 
tangular  patch  with  level  /—  1  grids  and  by  filling  the  overlapping  pieces.  However, 
it  may  be  necessary  to  go  to  even  coarser  grids  to  supply  these  level  I  -  1  values. 
This  is  done  by  applying  (i)  to  (iii)  recursively  to  the  smallest  rectangular  patch 
containing  the  unfilled  cells. 

For  efficiency,  it  is  important  that  the  boundary  values  are  supplied  on  a  rec¬ 
tangular  patch  at  one  time  and  not  computed  a  point  at  a  time.  Even  though  the 
amount  of  work  is  proportional  to  the  boundary  of  each  grid,  our  initial  implemen¬ 
tation  took  40  °o  of  the  run  time  and  had  to  be  rewritten.  By  working  on  grid 
patches,  the  bulk  of  the  memory  transfers  are  done  in  blocks,  and  the  number  of 
subroutine  calls  is  minimized.  This  is  particularly  important  on  the  Cray,  where 
there  is  a  substantial  performance  penalty  for  single  word  accesses  and  subroutine 
calls. 


5.  Creating  the  Grid  Hierarchy 

At  specified  time  intervals,  an  error  estimation  procedure  is  invoked,  and  a  new 
grid  structure  is  determined.  If  there  are  several  nested  levels  of  refined  grids,  the 
error  estimation  and  grid  generation  procedures  are  recursively  applied  on  each 
level,  from  finest  to  coarsest,  to  (re-)create  the  fine  grids  at  the  next  level.  The  error 
estimation  procedure  (see  Section  6)  produces  a  list  of  coarse  grid  points  with  large 
error  estimates,  indicating  that  a  fine  grid  patch  is  needed  in  that  region.  Every 
flagged  coarse  grid  point  should  be  included  in  a  finer  grid.  Our  grid  generation 
algorithms  try  to  produce  grids  that  have  as  little  overlap  as  possible,  so  that  the 
area  that  is  unnecessarily  refined  is  as  small  as  possible.  The  algorithm  also  strives 
for  a  small  number  of  patches  that  are  as  large  as  possible,  to  reduce  the  computa¬ 
tional  overhead.  It  is  difficult  to  find  a  foolproof  algorithm  that  satisfies  these  often 
conflicting  goals.  However,  we  have  developed  heuristics  that  have  been  successfully 
tested  in  many  different  applications.  A  much  fuller  discussion  of  grid  generation  is 
in  [3],  Here,  we  will  describe  the  particular  set  of  algorithms  that  produced  the 
numerical  results  in  Section  7. 

Suppose  there  is  a  base  level.  /ba(e,  where  grids  will  stay  fixed,  but  that  the  finer 
levels  from  lhi,e+  1  to  /rme„  may  be  "moved."  Starting  with  the  finest  level  grids,  we 
estimate  the  error,  using  a  procedure  described  next.  If  there  are  points  where  the 
error  estimate  is  too  high,  these  points  arc  flagged,  and  a  level  +  I  grid  will  be 
needed.  Next,  we  estimate  the  error  on  the  existing  -  t  grids.  If  there  are 
flagged  points,  a  different  level  /nnesl  grid  will  be  created,  making  sure  that  if  there 
are  any  level  /rinc„  +  1  grids,  they  are  properly  contained  in  the  level  /„„„,  grids.  This 
continues  until  the  error  is  estimated  on  the  base  level  grids.  Thus,  it  is  only 
possible  to  add  one  new  level  at  a  time,  although  many  levels  may  be  removed 
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during  a  single  regridding  operation.  (At  the  initial  time  however,  where  the  initial 
data  is  known  for  all  v,  v  and  not  just  at  coarse  grid  points,  it  is  possible  to  add 
many  levels  at  a  time.  This  is  essential  for  some  problems,  where  the  error  can 
depend  entirely  on  the  initial  conditions.  I 

In  more  detail,  our  regridding  algorithm  performs  the  following  steps: 

1 1 )  Adds  the  buffer  zone.  A  buffer  zone  of  unflagged  points  is  added  around 
every  grid.  This  ensures  that  discontinuities  or  other  regions  of  high  error  do  not 
propagate  out  from  a  fine  grid  into  coarser  regions  before  the  next  regridding  time. 
This  is  possible  because  of  the  finite  propagation  speed  of  hyperbolic  systems.  The 
larger  the  buffer  zone,  the  more  expensive  it  is  to  integrate  the  solution  on  the  fine 
grids,  but  the  less  often  the  error  needs  to  be  estimated  on  the  coarse  grids  and  the 
fine  grids  moved.  The  buffer  zone  is  added  by  flagging  all  coarse  grid  points  that 
are  sufficiently  close  to  flagged  points  with  high  error  estimates.  A  buffer  zone  of 
two  cells  in  each  direction  is  typical.  By  flagging  neighboring  points,  instead  of 
enlarging  grids  at  a  later  step,  the  area  of  overlap  between  grids  is  reduced. 

1 2  I  Flags  every  cell  at  level  l  corresponding  to  an  interior  cell  in  a  level  !  -  2 
grid.  This  will  maintain  proper  nesting,  by  ensuring  that  there  will  be  a  new  level 
/  -  1  grid  containing  every  point  in  the  level  /  *  2  grid,  even  if  the  level  /  grid  error 
estimation  did  not  report  a  high  error.  This  procedure  ensures  that  the  fine  grid 
error  estimates  are  used  instead  of  the  coarse  grid  estimates  at  the  same  point.  To 
ensure  proper  nesting,  points  within  one  cell  of  a  non-physical  (interior)  boundary 
of  G  are  deleted  from  the  list  of  flagged  points. 

(3 1  Creates  rectangular  fine  grids.  The  grid  generator  takes  all  the  flagged 
points  as  input,  and  outputs  a  list  of  corners  of  rectangles  that  are  the  level  /  -  1 
grids.  Nearby  points  are  clustered  together,  and  a  fine  grid  patch  spanning  each 
cluster  is  formed.  These  clustering  algorithms  use  heuristic  procedures  described 
separately  below. 

(4|  Ensures  proper  nesting.  The  new  fine  grids  are  checked  to  ensure  that 
they  are  properly  contained  in  the  base  level  grids.  If  they  are  not.  the  new  grid  is 
repeatedly  subdivided  until  each  piece  does  fit.  Since  the  flagged  points  originally 
were  inside  the  base  grid,  at  least  one  cell  from  the  boundary,  the  new  grid  contain¬ 
ing  the  flagged  points  must  eventually  lie  inside  as  well  Since  the  base  level  grids 
did  not  move,  step  (21  cannot  be  used  to  ensure  the  proper  nesting  of  this  ievef 
This  problem  only  arises  when  the  base  grids  are  a  non-convex  union  of  rectangles. 

Step  (3)  is  the  difficult  one.  Since  problems  in  gas  dynamics  develop 
1-dimensional  discontinuities,  we  have  streamlined  the  more  general  grid  generation 
procedures  of  [3]  for  this  particular  application.  The  procedure  we  use  here 
includes  a  bisection  step  and  a  merging  step.  Initially,  a  grid  patch  is  formed  around 
the  entire  list  of  flagged  cells  on  a  given  level.  The  efficiency  of  the  patch  is 
measured  by  taking  the  ratio  of  flagged  cells  to  the  total  number  of  cells  in  the  new 
grid.  If  this  efficiency  rating  is  less  than  an  input  minimum  efficiency  leg..  6U°oi. 
the  long  direction  of  the  rectangular  grid  is  bisected,  and  the  flagged  points  are 
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sorted  into  two  clusters  depending  on  which  half  they  are  in.  The  process  is 
repeated  on  the  two  clusters.  The  bisection  steps  ends  when  each  cluster  has  an 
acceptable  efficiency  rating. 

The  bisection  step  uses  no  geometric  information,  so  although  each  grid  may  be 
"acceptable"  by  itself,  the  resulting  grid  hierarchy  may  not  be  optimal.  For  this 
reason,  the  bisection  is  followed  by  a  merge  step.  In  addition  to  an  absolute 
efficiency  criterion,  grids  are  merged  if  the  new  grid  is  relatively  more  efficient  than 
the  two  smaller  grids.  The  cost  function  we  use  to  measure  this  is  proportional  to 
the  cost  of  an  integration  step  on  each  grid.  On  an  m  by  n  grid,  (m  •+•  1 )  by  <«  ■+■  1 1 
fluxes  are  calculated,  with  perhaps  1000  vector  operations  per  flux.  In  addition 
there  is  a  cost  associated  with  the  perimeter  of  each  grid:  finding  interface  condi¬ 
tions.  conservative  updating  of  coarser  grids,  and  special  slope  calculations  that  are 
done  only  for  boundary  fluxes.  Some  of  this  work  uses  scalar  arithmetic,  at  least  on 
machines  such  as  the  Cray  1  that  does  not  vectorize  indirect  addressing  and 
gather  scatter  operations.  The  total  cost  associated  with  a  grid  is  proportional  to 
»m  +  m  -r  n.  Grids  are  merged  if  the  single  resulting  grid  has  a  smaller  cost.  The 
merging  step  ends  when  no  pair  of  grids  can  be  successfully  merged. 

Although  this  procedure  is  somewhat  ad  hoc,  it  has  been  successfully  used  on 
several  different  types  of  problems.  The  grid  generation  routines,  not  including  the 
solution  initialization  on  each  grid  or  the  error  estimation  to  produce  the  flagged 
points,  account  for  approximately  1.7 °b  of  the  CPU  time  for  a  typical  run. 


6.  Error  Estimation 

In  [6],  estimates  of  the  local  truncation  error  were  used  to  select  those  grid 
points  on  a  given  level  with  unacceptably  large  errors.  If  the  solution  u(.v.  n  is 
smooth  enough,  the  local  truncation  error  u( x.  i  +  k\  —  Qu{ x.  1 1  on  a  mesh  with 
spatial  step  h  and  time  step  k  satisfies 

ui  a.  r  ~  k  i  -  Qu(  v.  n  =  Ar(c,(.v.  t \k*  -t-  cax.  t  l/C)  ■+■  kO(k v  * 1  +  h*  * 1 ) 
srl.v.  t  la-  kOik*  ''+/>■''  1 ). 

where  the  leading  term  is  denoted  by  r.  Here  we  assume  our  difference  method  Q 
has  order  of  accuracy  q  in  both  time  and  space.  If  u  is  smooth  enough,  then  if  we 
take  two  time  steps  with  the  method  Q.  to  leading  order  the  error  is  2r  . 

in  v.  i  ~  2k)  -  Q-m x.  i)  =  2r  ■+■  kO{ k*  *  :  +•  hu ' 1 ). 

Let  Q:i,  denote  the  same  difference  method  as  Q  but  based  on  a  mesh  widths  of  2 h 
and  2k.  Then 


ui  v.  r  —  2A-  >  —  Q:i,  ul  v.  rl  =  2*  *  'r  -*-  Oi/C  * :  |. 
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By  taking  two  steps  with  the  regular  integration  scheme,  and  one  "giant''  step  usu 
every  other  grid  point,  the  difference 


Q'~u[  v.  1 1  -  Q:,.ut  x.  1 1 


=  r  -Oi/r  "-I 


gives  an  estimate  of  the  local  truncation  error  at  time  r.  We  emphasize  that 
necessary  to  know  the  exact  form  of  the  truncation  error  te.g..  Iru  <. 
order 

This  procedure  is  easily  implemented  in  AMR  for  the  conservative 
ference  methods  presented  here.  The  values  on  a  grid  at  a  given  level  are 
onto  a  virtual  grid  coarsened  by  a  factor  of  two  in  each  direction.  The  -.m 
both  grids  is  ad'.anced  in  time:  the  original  grid  for  two  time  sims.  the  .. 
grid  for  one  step  using  a  time  step  twice  as  large.  The  difference  between  ' 
tions  obtained  on  the  two  grids  at  each  point  is  proportional  to  the  loc.u  r 
error  at  that  point.  At  coarse  cells  where  the  difference  between  the  :.v 
values  exceed  some  tolerance,  all  four  cells  contained  in  the  real  grid  are 
requiring  refinement.  Notice  that  this  estimation  procedure  is  tndeponder 
finite  difference  method  actually  used,  as  well  as  the  pde.  One  disadv.int  ig 
procedure  is  that  it  always  predicts  a  large  error  in  the  neighborhood  . 
discontinuities,  ft  is  easy  to  construct  examples  for  wb-.a  the  procedv.re 
above  will  give  values  on  the  coarsened  grid  whwii  differ  pointwoe  bv 
independent  of  the  mesh  spacing  in  the  neighborhood  of  a  shock  In  gc:.. 
leads  to  refinement  of  the  mesh  at  all  discontinuities  with  strength  c. 
some  minimum. 

Theoretically,  one  could  define  a  distributional  error  by  averaging  m. 
between  the  two  solutions  over  vome  region  centered  at  the  given  -eii 
is  Oi !  i  relative  to  the  mesh  spacing.  We  have  w  sed  various  technic...- 
mg  out  such  a  procedure  all  of  which  are  equivalent  to  ignoring  the  po .- " 
estimate  in  the  neighborhood  of  those  discontinuities  which,  bv  -onto 
are  considered  adequatelv  resolved.  For  problems  in  shock  hydrodv 
shock  discontinuities,  and  not  slip  surfaces  or  contact  discontinuities 
satisfv  anv  such  criteria.  This  is  because  conservative  finite  differ.!  .. 


applied  to  shocks  mimic  the  convergence  of  characteristics  in  the  anai.'  . 
so  that  the  shock  spreads  only  over  a  fixed  number  of  zones  mdeper.c.-  •  ..■ 
mesh  spacing  and  time.  Thus,  for  example,  it  is  unnecessary  to  reline 
neighborhood  of  a  shock  separating  two  constant  states  In  conira-:  '  . 

necessary  to  refine  at  linearly  degenerate  discontinuities  since  the  tuirrmc-  .  .i- 
over  which  they  spread  is  an  increasing  function  of  time. 

There  is  a  second  set  of  difficulties  with  refinement  in  the  preserve  •  -g 
shocks  There  is  evidence  that  shock-capturing  methods  are  zeroth-order  S. 

i  e  .  ih.tt  the  fluxes  computed  in  the  neighborhood  of  the  shock  differ  m  <  >  ■  m 

the  exact  fluxes  at  a  given  time  step  [21.  16],  These  Oil  i  errors  could  j.-  ■ 
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w-^ves  associated  with  the  characteristic  families  crossing  the  shock,  sending  Ot  I ) 
pointwise  errors  into  the  postshock  region.  For  shocks  computed  on  a  single 
uniform  grid,  this  does  not  occur,  because  the  same  0(  1 )  errors  are  committed  on 
successive  cell  edges  with  a  phase  lag,  so  that  errors  in  the  time  integral  of  suc¬ 
cessive  fluxes  cancel  upon  differencing.  However,  when  a  shock  intersects  an  inter¬ 
face  between  two  grids  with  different  mesh  spacings.  the  0(  1  |  errors  in  the  time 
integrals  of  the  fluxes  on  each  of  the  two  grids  will  be  different,  generating  ()\  \  i 
errors  in  the  solution  propagating  into  the  postshock  state  In  practice,  we  have 
observed  that  the  amplitude  of  the  spurious  waves  generated  in  this  fashion  is 
proportional  to  the  amplitude  of  the  jump  in  the  characteristic  quantities  carried  by 
characteristics  crossing  the  shock.  In  practice,  then,  there  is  usually  a  threshold 
shock  strength  below  which  the  errors  generated  by  a  shock  crossing  a  grid  discon¬ 
tinuity  are  acceptable  and  above  which  the  errors  generated  are  too  large.  For  the 
latter  shocks,  they  must  be  refined  everywhere,  if  they  are  to  be  refined  unyw  iere 

We  illustrate  this  with  an  example  where  we  force  the  algorithm  not  to  refine  the 
grid  abi  ve  a  certain  height.  This  forces  the  strong  incident  shock,  with  a  shock 
Mach  number  of  10.  to  pass  through  a  fine  grid  boundary  into  the  coarse  grid.  The 
oscillations  caused  by  this  are  apparent  in  the  contour  plots  of  Fig.  6.1  and  the  plot 
of  Fig.  6.2  for  a  fixed  value  of  i 

Combining  the  two  considerations  given  above,  a  fairly  general  refinement 
strategy  is  to  use  the  local  truncation  error  estimate  described  above,  but  ignore  it 
in  the  neighborhood  of  gus-dynamtc  shocks  whose  strength  lies  below  some 
predetermined  threshold.  This  has  the  effect  of  refining,  possibly  unnecessarily.  all 
shocks  whose  strength  is  above  the  threshold.  We  have  found  that  this  strategy 
works  acceptably  well  if  the  coarsened  base  grid  is  sufficiently  fine,  so  that  the 
waves  are  well  separated.  A  much  simpler  strategy,  which  is  applicable  in  a  large 
class  of  problems,  is  simply  to  use  the  user's  knowledge  of  the  problem  instead  oi 
the  truncation  error  estimates.  For  example,  in  the  shock  reflection  problems  giver, 
below,  the  solution  is  made  up  entirely  of  smooth  waves  and  weak  shocks  a  certain 
distance  behind  the  incident  shock.  Consequently,  we  simply  ignore  the  local 
truncation  error  estimate  in  that  region. 


Fig  f)  !  Contour  plot  showing  the  effects  of  a  strong  shock  passing  through  a  end  boundarv 
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Fit i  ft  2  Density  profile  for  a  horizontal  line  slightly  below  the  grid  interface  intersecting  the  -n.^k 


~  Ni  merical  Results 

We  choose  as  our  test  problem  reflection  of  a  planar  shock  in  an  ideal  gas  h>  an 
oblique  surface.  In  this  problem,  a  straight  shock  is  incident  on  a  perfectly  reflecting 
surface.  At  later  times,  a  reflected  wave  pattern  is  generated,  depending  only  on  2, 
A/,.  and  ,  .  where  a  is  the  angle  between  the  direction  of  propagation  and  the  reflec¬ 
ting  surface.  A/,  is  the  incident  shock  Mach  number,  and  ;  is  the  ratio  of  specific 
heats.  The  solution  to  this  problem  is  formally  self-similar,  depending  on  n.  ..  n  1 
only  in  the  combination  u  r.  1  n.  Thus,  the  time  dependence  of  the  solution  :>  1 
given  by  a  linear  scaling  of  a  tlxed  wave  pattern  with  time.  We  are  particularly  i 
interested  in  values  of  the  problem  parameters  for  which  very  complicated  small  ‘ 
scale  structures  are  observed. 

The  underlying  integration  method  in  our  AMR  calculations  is  a  second-order 
Godunov  method  described  in  [10].  In  our  calculations,  we  make  use  of  the  fact 
that  the  regions  where  the  small  scale  behavior  can  appear  are  localized  in  the 
vicinity  of  the  reflection  point.  Our  procedure  for  estimating  the  error  is  to  measure 
the  local  truncation  error  of  the  density,  except  that  we  do  not  tag  points  beyond 
a  certain  distance  behind  the  incident  shock.  This  effectnelv  restricts  our  refinement 
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TABLE  I 


Breakdown  of  Computational  Times  Obtained  Using  FLOWTRACE 


Grid  integration 

78.29 

Interpolation 

12.95 

Output 

2.87 

Grid  updates 

2.’8 

Grid  generation 

i.7! 

Memory  management 

0.59 

to  be  in  a  window  moving  with  the  reflection  point,  and  shuts  off  the  grid  refine¬ 
ment  ground  the  iweuki  reflected  shock. 

The  results  presented  here  were  calculated  on  the  Cray  XMP  22  at  the  LLNL 
NMFECC.  using  the  CFT  1.14  compiler.  To  obtain  detailed  diagnostic  information 
about  where  most  of  the  time  is  spent  in  the  calculation,  we  used  the  FLOW- 
TRACE  option  of  the  CFT  compiler  in  the  first  calculation  below.  The  total  time 
spent  was  5674  seconds  of  CPU  time.  Table  I  shows  a  breakdown  of  the  calculation 
time  into  six  categories:  the  integration  routine,  the  interpolation  routines  (for  con¬ 
structing  boundary  conditions  and  initializing  new  fine  grids),  the  updating  routines 
(fine  grids  updating  coarse  grids  and  for  maintaining  conservation  across  grid 
interfaces),  the  grid  generation  routines,  output  routines,  and  memory  management 
routines. 

The  main  result  is  that  the  integration  step  takes  about  80%  of  the  computa¬ 
tional  time.  This  figure  includes  integration  steps  needed  for  the  error  estimation. 
However,  measurements  show  that  the  latter  is  only  3%  of  the  integration  cost, 
with  actual  useful  integration  steps  accounting  for  97%  of  the  integration  time. 
Note  that  the  error  estimation  cost  is  very  small  despite  the  fact  the  error  is 
estimated  at  every  other  coarse  time  step.  There  are  two  reasons  for  this.  First,  over 
90 "o  of  the  cells  being  integrated  belong  to  the  finest  level  grids  (level  3  in  both 
calculations  i.  and  the  error  is  not  estimated  there.  Second,  since  refinement  is  per¬ 
formed  in  time  as  well  as  space,  the  overwhelming  majority  of  the  work  is  done  on 
the  finest  grid.  Table  II  shows  the  number  of  cell  updates  done  on  each  grid  level, 
as  well  as  the  total  number  of  cell  updates  done  for  error  estimation. 

We  can  thus  obtain  a  rough  estimate  of  the  efficiency  of  AMR  relative  to  com¬ 
puting  on  a  uniform  grid.  About  80  %  of  the  run  time  is  spent  integrating  grids  The 

TABLE  II 


Number  of  Cell  Updates  at  Each  Level 


Level  I 

1  J 

v£> 

oo 

O 

Level  2 

4  59  *  10" 

Level  3 

1  13  *  10' 

Error  estimation 

3  06  «  10" 
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finest  level  grids  occupy  only  about  10%  of  the  domain.  Thus  an  equivalent 
uniform  grid  computation  would  require  a  factor  of  8  more  CPU  time.  Of  course 
in  this  particular  problem,  one  could  omit  computations  ahead  of  the  incident 
shock,  since  the  solution  there  is  constant,  saving  approximately  half  the  uniform 
grid  time.  In  general,  this  could  not  be  done. 

It  is  more  difficult  to  compare  memory  usage  with  that  of  a  uniform  grid  calcula¬ 
tion.  The  integration  algorithm  used  here  requires  five  2-dimensional  grid  arrays  for 
scratch  space,  in  addition  to  the  four  required  to  store  the  conserved  quantities, 
so  that  the  memory  requirements  for  a  uniform  grid  calculation  would  be 
9  ■  320  1600  ^  4.5  x  106  words.  In  contrast,  the  maximum  storage  used  in  the 
calculation  performed  here  was  8.94  x  10?  words.  Much  of  the  memory  use  in  -\Nf  R 
is  due  to  saving  two  time  levels  of  the  solution  on  each  grid.  It  is  possible  to  avoid 
the  memory  overhead  of  having  full  grid  scratch  arrays  by  breaking  the  calculation 
into  pieces.  (In  fact,  we  effectively  do  this  in  the  AMR  calculations  by  restricting  the 
size  of  any  grid  to  be  less  than  some  pre-determined  maximum,  subdividing  grids 
that  are  too  large.)  However,  this  would  introduce  overheads  and  programming 
complexities  in  the  uniform  grid  calculation  similar  to  those  in  AMR  In  any  case, 
even  if  those  overheads  could  be  neglected  and  only  four  full  grid  arrays  were 
required,  the  memory  required  would  be  2.0  x  105  words,  a  factor  of  2.2  larger  than 
that  required  by  AMR. 

In  Fig.  7.1,  we  show  results  for  the  case  V/s  =  10.  x  =  30  ,  ;  =  1.4.  The  domain  is 
a  rectangle  of  length  2.0  by  0.4.  with  initial  coarse  grid  spacing  J.v  =  J 1  =  0 1 P  The 
calculation  ran  for  149  coarse  grid  time  steps.  The  error  was  estimated  every  other 
step,  with  a  buffer  zone  of  one  cell  and  a  grid  efficiency  of  65  %.  The  error  tolerance 
was  0.02.  The  mesh  was  refined  by  a  factor  of  4  in  each  direction  at  each  grid  level 
The  finest  grids  in  this  calculation  represent  a  factor  of  4  increase  in  resolution  in 
each  spatial  direction  over  the  finest  grid  calculation  in  [22].  Figure  "  la  show*  'he 
location  of  the  level  2  and  3  grids  at  time  r=  1.20.  In  displaying  the  solution,  .ve 
show  two  sets  of  plots.  One  is  a  contour  plot  of  the  full  flow  field.  The  other  o  an 
enlargement  of  the  region  around  the  reflection  point.  This  is  the  part  ot  ; he 
domain  covered  by  the  level  two  and  three  grids.  In  both  cases,  the  contour  plots 
are  made  using  the  finest  available  grid  in  the  subregion.  Due  to  the  increu'ed 
resolution,  we  can  now  observe  a  non-self-similar  Kelvin-Helmholtz  rollup  aUrj 
the  principal  slip  line.  This  is  to  be  expected,  since  this  slip  line  is  instuhie  TLe 
Kelvin-Helmholtz  rolls  are  formed  near  where  the  weak  shock  emanating  from  tie 
second  triple  point  impinges  on  the  slip  line  They  then  propagate  along  the  v,:p 
line  and  are  eventually  swept  up  into  a  large  rollup  at  the  tip  of  the  jet.  along  "e 
bottom  wall. 


Finally,  in  Fig.  7.2  we  present  results  for  \fs  =  8.  x  =  35  .  ;  =  1.107  It  h.iv  ^eea 
noticed  [11]  that  the  wave  patterns  associated  with  double  Mach  reflection  become 
increasingly  complex  as  y  is  reduced.  The  jet  along  the  reflecting  wall  formed  by  ;h 


slip  line  from  the  principal  Mach  triple  point  is  more  and  more  •orongi* 
accelerated,  pushing  the  Mach  stem  out  ahead  of  it.  This  leads  to  strongly  rota- 
tional  supersonic  flow  and  the  formation  of  multiple  Mach  triple  points  The  pre- 


LOCAL  ADAPTIVE  MESH  REFINEMENT 


Fig.  I —  Continued 


LOCAL  ADAPTIVE  MESH  REFINEMENT 


83 


sent  results  represent  a  rather  extreme  example  of  this,  with  a  total  of  seven  Mach 
triple  points  in  the  double  Mach  region,  including  a  third  triple  point  along  the  top 
shock.  The  seven  triple  points  are  marked  in  Fig.  7. 2d.  This  calculation  gives  some 
indication  of  why  adaptive  mesh  refinement  is  an  important  tool  in  these  problems. 
As  7  approaches  one.  the  distance  between  the  leading  edge  of  the  wall  jet  and  the 
main  Mach  stem  becomes  smaller  and  smaller,  requiring  more  grid  resolution  in 
that  region.  The  value  of  y  in  this  calculation  represents  the  limit  of  resolution, 
given  the  resources  available  on  the  Cray  XMP.  The  calculations  in  [II]  using 
uniform  grids  without  mesh  refinement,  were  not  fully  resolved  for  ;<  1.25.  Even 
with  mesh  refinement,  at  this  low  value  of  y  the  flow  field  is  not  fully  resolved  at 
time  /  =  0.115.  in  Fig.  7.2a  and  b.  Only  after  running  to  time  t  =  0.230.  which  by 
self-similarity  corresponds  to  increased  grid  resolution,  is  the  solution  adequately 
resolved. 


8.  Conclusions 

The  complexity  of  our  AMR  code  might  be  intimidating  to  a  new  user.  Not 
counting  the  integration  routine,  our  program  consists  of  3000  lines  of  Fortran. 
However,  a  big  code  is  not  necessarily  a  fragile  code.  We  have  been  careful  to 
develop  AMR  to  make  it  automatic  and  robust.  In  addition,  a  user  should  be  able 
to  use  AMR  without  having  to  understand  it  all.  This  makes  it  important  to 
develop  AMR  in  a  modular  way.  A  user  should  be  able  to  plug  in  an  integrator  for 
a  new  problem  without  knowing  details  about  how  the  more  computer  science 
oriented  parts  of  AMR  work,  but  knowing  that  these  other  parts  will  indeed  work. 
We  have  already  demonstrated  this  modularity  by  using  AMR  to  compute  tran¬ 
sonic  flow  in  conjunction  with  FL052  [5]  and  to  compute  a  combustion  problem 
with  a  simple  induction  time  model  for  chemistry  [2]. 

The  most  difficult  problems  will  best  be  solved  by  combining  several  adaptive 
techniques.  Despite  its  more  complicated  data  structures.  AMR  has  already  been 
combined  with  the  conservative  front-tracking  scheme  of  [9],  This  enables  tracking 
of  a  strong  incident  shock,  while  using  shock-capturing  for  the  other  discontinuities, 
and  avoids  the  mismatch  of  strong  captured  shocks  crossing  grid  boundaries. 
This  combined  approach  is  being  used  to  study  transition  from  regular  to  Mach 
reflection.  Finally,  we  intend  to  couple  this  method  with  the  variational  technique 
of  [8].  Their  mesh-moving  technique  would  allow  the  underlying  mesh  geometry  to 
be  approximately  aligned  with  global  features  in  the  flow,  leading  to  more  efficient 
refined  meshes.  However,  the  actual  mesh  refinement  for  error  reduction  would  be 
done  with  AMR.  so  the  global  time  step  penalty  of  moving  mesh  methods  is  not 
incurred.  Lastly,  a  major  open  question  is  how  to  use  implicit  difference  schemes 
with  embedded  grids  for  a  time-dependent  calculation.  This  will  be  needed  to 
compute  solutions  to  hyperbolic-parabolic  problems,  such  as  the  Navier-Stokes 
equations  at  high  Reynolds  number. 
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Abstract 

We  present  a  Cartesian  mesh  algorithm  with  adap¬ 
tive  refinement  to  compute  flows  around  arbitrary 
geometries.  Cartesian  meshes  have  been  less  popular 
than  unstructured  or  body-fitted  meshes  because  of 
several  technical  difficulties.  We  present  an  approach 
that  resolves  many  of  these  problems.  Cartesian  meshes 
have  the  advantage  of  allowing  the  use  of  high  resolution 
methods  that  are  difficult  to  develop  on  unstructured 
grids.  They  also  allow  for  efficient  implementation  on 
vector  computers  without  using  gather-scatter  operations 
'except  at  boundary  cells.  Some  preliminary  computa¬ 
tional  results  using  lower  order  boundary  conditions  are 
presented. 

1.  Introduction 

The  construction  of  logically  rectangular  body- 
fitted  grids  for  complicated  geometries  is  notoriously 
difficult.  One  alternative  is  to  use  an  unstructured  mesh, 
so  that  the  cell  volumes  are  not  derived  by  a  smooth 
mapping  from  a  rectangular  domain,  as  in  [16]  for  exam¬ 
ple.  Another  possibility  is  to  simply  use  a  Cartesian  grid 
over  the  entire  flowfield.  This  introduces  the  difficulty 
of  imposing  solid  wall  boundary  conditions  on  a  grid  that 
is  not  aligned  with  the  body. 

Nonetheless,  there  are  several  reasons  to  prefer  the 
Cartesian  grid  approach,  in  addition  to  the  ease  of  grid 
generation.  First,  it  allows  the  use  of  higher  order  accu¬ 
rate  shock  capturing  methods  that  are  difficult  to  achieve 
on  an  unstructured  mesh  with  no  coordinate  directions. 
A  Cartesian  grid  integrator  is  highly  vectorizable. 
Gather-scatter  type  vector  operations  need  to  be  per¬ 
formed  only  in  a  lower  dimensional  region,  and  not  over 
the  entire  flowfield.  The  basic  solver  on  a  Cartesian 
mesh  is  also  simpler  than  on  a  body-fitted  mesh,  since 
there  are  no  metric  terms.  Finally,  there  is  some  evi¬ 


dence  [20]  that  for  strong  shock  calculations  an  unstruc¬ 
tured  mesh  has  larger  phase  errors,  and  thus  poorer 
shock-capturing  abilities,  than  structured  grids. 

There  are  also  several  difficulties  in  using  Carte¬ 
sian  grids.  The  main  difficulty  is  the  small  cell  problem. 
Arbitrarily  small  cells  arise  at  the  edge  of  the  domain 
where  the  grid  intersects  a  body.  Stable,  accurate  and 
conservative  difference  schemes  are  needed  for  these 
cells.  We  would  like  the  time  step  for  a  time-accurate 
computation  to  be  based  on  the  cell  volume  of  the  regu¬ 
lar  cells  away  from  the  body,  and  not  be  restricted  by 
small  cells  at  the  boundary.  The  time  step  appropriate  for 
the  regular  grid  cells  can  give  a  Courant  number  that  is 
orders  of  magnitude  larger  than  1  for  the  smaller,  irregu¬ 
lar  cells.  This  will  lead  to  stability  problems  with  stan¬ 
dard  explicit  methods.  In  this  work,  we  use  an  approach 
based  on  wave  propagation  that  essentially  increases  the 
size  of  the  stencil  near  these  small  cells  and  maintains 
stability  for  arbitrary  time  steps.  This  large  time  step 
approach  was  studied  in  one  space  dimension  in  [13], 
and  has  been  applied  in  the  present  context  of  small 
boundary  cells  in  [14,15]. 

Another  problem  with  Cartesian  grids  is  one  of 
accuracy.  Grid  stretching,  used  in  body-fitted  grids  to 
cluster  the  grid  points  in  regions  where  they  arc  needed, 
cannot  be  used.  Moreover,  since  the  grid  is  not  aligned 
with  the  boundary,  a  loss  of  resolution  may  occur  near 
the  boundary.  To  improve  the  accuracy  we  use  an  adap¬ 
tive  mesh  refinement  algorithm  developed  in  [8].  Rec¬ 
tangular  refined  grids  are  superimposed  on  the  coarse 
grid,  so  that  the  efficiency  of  the  integrator  on  each  grid 
is  maintained.  The  time  step  is  refined  along  with  the 
mesh  width  on  the  fine  grids,  so  that  the  CFL  condition  is 
maintained  while  allowing  larger  time  steps  on  the 
coarser  grids.  This  further  concentrates  the  computa¬ 
tional  work  where  it  is  needed. 
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Cartesian  mesh  methods  have  received  increased 
attention  recently.  Cartesian  grids  were  used  in  [19]  to 
solve  the  full  potential  equations.  This  was  extended  to 
the  Euler  equations  in  two  space  dimensions  in  [1 1],  and 
to  three  dimensions  in  [12].  Cartesian  grids  have  also 
been  used  in  conjunction  with  an  implicit,  flux-vector 
split  method  for  the  Euler  equations  [10].  These  calcula¬ 
tions,  however,  suffer  from  the  lack  of  resolution  of  a 
Cartesian  mesh;  none  of  the  calculations  used  a  local 
mesh  refinement  algorithm.  The  use  of  a  global,  tensor 
product  grid  to  concentrate  grid  lines  near  the  leading 
edge  of  an  airfoil,  for  example,  can  be  very  wasteful.  In 
those  computations,  the  small,  irregular  cells  near  the 
body  were  merged  into  their  neighboring  cells  to  create  a 
cell  that  was  large  enough  to  satisfy  a  stability  constraint. 
This  procedure  loses  resolution. 

In  the  next  section,  we  describe  our  overall  com¬ 
putational  method,  including  the  organization  of  the  cal¬ 
culation  and  the  boundary  representation.  Section  3 
describes  the  large  time  step  method  and  the  implemen¬ 
tation  of  the  solid  wall  boundary  conditions.  Section  4 
describes  the  mesh  refinement  algorithm.  In  Section  S 
we  show  applications  of  our  method  to  the  inviscid  Euler 
equations  in  two  geometries;  shocked  flow  around  two 
cylinders,  and  a  curved  channel  calculation. 

2.  Algorithm  and  Data  Structures 

We  consider  the  Euler  equations  in  two  space 
dirfiensions, 

U,  +  F(U)t  +  G(U)y  =  0, 
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P  =  (Y-l)(p£-p(u2+v2)/2)). 

The  boundary  of  a  solid  object  is  approximated  by 
piecewise  linear  segments  at  shown  in  Figure  la.  We 
assume  that  the  boundary  cuts  through  each  cell  at  most 
once,  so  that  the  resulting  irregular  cell  is  a  polygon  with 
at  most  five  sides.  We  use  a  finite  volume  method,  with 
Ujj  representing  the  cell  average  of  the  vector  of  con¬ 
served  quantities  in  cell  i.e. 

where  Ay  is  the  area  of  the  cell.  For  regular  cells,  away 
from  the  boundary.  Ay  =  Ax  Ay,  but  Ay  may  be  orders  of 
magnitude  smaller  for  cells  on  the  boundary. 


All  grid  cells  are  indexed  using  the  rectangular 
Cartesian  structure.  Additional  information  about  the 
irregular  cells  is  kept  in  a  linked  list  data  structure  that  is 
easily  traversed  in  implementing  the  boundary  condi¬ 
tions.  Necessary  information  includes  the  cell  area  and 
list  of  vertices  for  each  irregular  cell,  as  well  as  a  pointer 
to  its  location  in  the  Cartesian  grid.  In  the  other  direc¬ 
tion,  a  two  dimensional  integer  array  indicates  whether  a 
Cartesian  cell  is  regular,  and  for  irregular  cells  contains  a 
pointer  to  the  corresponding  location  in  the  linked  list. 
When  grid  refinement  is  used  there  may  be  several  rec¬ 
tangular  grids,  each  with  its  own  solution  storage  and 
irregular  points  list. 


(a) 


^d+1/2 


(b) 

Figure  1  (a)  a  Cartesian  grid.  The  shaded  region 
represents  a  solid  body,  (b)  Blowup  of  cell  (i.y) 
showing  the  fluxes. 

In  each  time  step,  the  cell  values  are  updated  by 
differencing  fluxes  at  the  cell  sides,  as  illustrated  in  Fig¬ 
ure  lb.  The  updating  formula  is 

Utf  =^-fL  [FMi  -  F^.y  (1) 

+  G«./4l/2  -  Gy-\u  +  Hy], 

boundary  of  the  cell.  For  example,  with  Godunov’s 


2 


method  we  would  obtain  Fi+1/2,y  by  solving  a  Riemann 
problem  u,  +  f  (u)x  =  0  with  left  and  right  states  Uy  and 
Ui+ 1  j  to  determine  the  correct  intermediate  state  u  .  We 
then  set  Fi+VXj  =  hi+VXjf(u*),  where  hi+u2J  is  the 
length  of  the  interface  between  cells  (i,y)  and  (i+lj). 
For  regular  cells  this  length  is  just  Ay.  Similarly,  G;j+  m 
is  the  flux  per  unit  time  through  the  top  of  the  cell. 
Finally,  Ht  J  is  the  flux  per  unit  time  through  the  irregular 
side  of  the  cell,  which  represents  the  solid  wall  boundary 
of  the  fluid  domain.  In  regular  cells,  Htj  =  0.  With 
higher  order  Godunov  schemes  such  as  MUSCL,  the  left 
and  right  states  are  modified  using  slope  information  to 
achieve  second  order  accuracy. 

In  irregular  cells  there  are  two  difficulties  with  this 
approach.  First,  the  neighboring  cell  values  needed  to 
define  appropriate  slopes  may  not  be  present.  This  is 
currently  handled  by  setting  the  slopes  to  zero,  so  that 
the  flux  reduces  to  the  Godunov  flux  at  these  interfaces. 
Improvements  to  this  algorithm  are  currently  under 
study. 

Even  with  first  order  fluxes,  there  is  still  a  stability 
problem.  Use  of  these  fluxes  in  updating  formula  (1) 
will  give  instabilities  in  cells  where  Atj  is  very  small 
relative  to  A/.  A  wave  propagation  interpretation  of  this 
is  given  in  Section  3,  where  we  present  a  way  to  modify 
the  fluxes  to  account  for  the  reflection  of  waves  at  the 
boundary  and  obtain  a  much  more  stable  algorithm. 
First  we  present  an  outline  of  the  overall  algorithm. 

Step  1 :  Initial  flux  computation.  In  the  first  pass, 
fluxes  at  all  cell  boundaries  are  computed  assuming  that 
the  grid  is  regular,  even  at  interfaces  where  both  neigh¬ 
boring  cells  lie  outside  of  the  actual  fluid  domain.  This 
is  done  for  ease  of  vectorization,  but  the  calculations  out¬ 
side  the  domain  will  have  no  influence  on  the  final  solu¬ 
tion.  Also,  the  interface  lengths  are  always  assumed  to 
be  Ax  or  Ay  in  computing  the  flux  per  unit  time,  regard¬ 
less  of  the  true  length  of  the  side.  This  will  be  corrected 
in  step  2  as  required. 

Step  2:  Flux  Modification  Near  the  Boundary.  In 
the  second  step,  we  march  around  the  solid  boundary  of 
the  fluid  domain,  modifying  the  fluxes  of  the  irregular 
cells.  First  we  adjust  the  fluxes  F  and  G  at  each  interface 
to  incorporate  the  correct  length  rather  than  the  standard 
length  Ax  or  Ay.  Next,  we  modify  the  fluxes  to  improve 
the  stability  of  the  small  cells  and  incorporate  the  solid 
wall  boundary  conditions.  This  will  be  described  in 
more  detail  in  the  next  section.  Finally,  we  calculate  the 
fluxes  at  the  irregular  side  of  each  boundary  cell. 
This  is  also  described  in  section  3. 

Step  3:  Updating  £/;,/.  The  grid  values  U,,y  are 
now  updated  using  the  flux  differencing  formula  (1). 

Step  4:  Smoothing  at  the  Boundary.  Although  the 
flux  modification  of  Step  2  is  intended  to  give  stability 


for  arbitrarily  small  cells,  in  practice  very  smalt  cells  still 
cause  difficulties  in  some  computations.  The  exact 
causes  of  this  are  currently  under  study.  In  the  present 
code,  stability  is  restored  by  means  of  an  averaging  pro¬ 
cess.  In  very  small  cells  (those  with  area  less  than  3%  of 
the  regular  cell  size,  typically),  the  value  of  U"*1  is 
replaced  by  a  weighted  average  of  the  original  value  and 
the  value  in  one  or  more  neighboring  cells.  The  choice  of 
cells  depends  on  the  local  geometry.  The  value  in  the 
neighboring  cetl(s)  is  also  modified  in  such  a  way  that 
conservation  is  maintained.  The  weights  used  arc  pro¬ 
portional  to  the  cell  areas  and  hence  the  value  in  the 
small  cell  is  replaced  by  a  value  that  is  essentially  equal 
to  that  of  the  cell’s  primary  neighbor  (or  a  weighted 
averaged  of  two  or  three  neighboring  cells).  This  pro¬ 
cedure  has  been  found  to  eliminate  any  remaining  insta¬ 
bilities.  This  algorithm  is  similar  to  the  flux  redistribution 
algorithm  in  [91, 

3.  Boundary  Conditions  and  Flux  Modification 

For  irregular  cells  at  the  boundary,  the  fluxes  that 
are  calculated  in  the  first  stage  of  the  algorithm  arc  sim¬ 
ply  the  Godunov  fluxes  obtained  by  solving  the  Riemann 
problem  between  this  cell  and  each  neighbor.  In  order  to 
understand  how  these  fluxes  should  be  modified  for 
small  cells,  we  first  consider  Godunov’s  method  on  a 
regular  grid  cell.  The  method  can  be  interpreted  in  the 
following  way.  Solve  the  Riemann  problem  at  each 
interface  to  obtain  waves  propagating  away  from  the 
interface.  For  each  wave  that  propagates  into  the  cell,  let 
A U  represent  the  jump  in  the  conserved  quantities  across 
the  wave.  Suppose  that  in  time  A t  this  wave  sweeps 
through  a  certain  fraction  a  of  the  cell.  Then  the  cell 
average  U,j  is  updated  by  the  quantity  a AU.  Note  that 
the  CFL  condition  requires  a  <  1,  and  that  a  is  the  ratio 
of  the  area  swept  out  by  the  wave  to  the  total  cell  area 
(see  Figure  2).  The  wave  shown  in  Figure  2  propagates 
with  speed  s  >0  from  the  left  side  of  the  cell,  and  so 

-  JA/Ay  _  sAl 
~  Nx  ‘ 

Now  consider  the  same  wave  but  suppose  that  the 
cell  in  question  is  an  irregular  cell  as  shown  in  Figure  3a. 
We  now  have  a  =  sNt  lAt  l  >  1,  and  updating  U^  by 
a  AU  would  lead  to  instability.  However,  an  alternative 
approach  that  is  physically  more  reasonable  is  to  update 
UtJ  by  only  l.OAU,  since  the  wave  overlaps  the  entire 
cell,  and  then  to  update  cells  further  to  the  right  by  the 
remainder  of  the  wave  (a-l)AU.  This  is  the  basis  of  the 
large  time  step  method  originally  described  in  [13].  In 
the  present  context  however,  there  are  no  cells  to  the 
right.  Instead,  there  is  a  solid  wall  boundary,  off  which 
the  wave  should  reflect.  This  is  illustrated  in  Figure  3b. 
The  portion  of  the  wave  that  lies  outside  the  domain  is 


3 


Ay 


Figure  2  The  wave  containing  the  jump  in  con¬ 
served  quantities  travels  to  the  right  a  distance  sAt 
away  from  the  interface. 


reflected  normal  to  the  boundary  segment  of  this  cell. 
Linearization  of  the  solid  wall  boundary  conditions  sug¬ 
gests  that  the  reflected  wave  should  carry  a  jump  A 0, 
which  has  the  same  jump  as  AU  in  density,  pressure  and 
tangential  velocity,  but  which  has  the  jump  in  normal 
velocity  negated.  (Normal  and  tangential  refer  to  the 
orientation  of  the  solid  wall  boundary). 

This  reflected  wave  overlaps  some  fraction  P<1  of 
the  cell,  and  so  Uu  is  further  updated  by  $A0.  In  addi¬ 
tion,  the  reflected  wave  may  overlap  neighboring  cells, 
and  each  of  these  is  also  updated  by  the  fraction  of  the 
cell  overlapped  multiplied  by  At/.  In  Figure  3b,  three 
neighboring  cells  are  affected  by  the  reflected  wave.  This 
is  the  maximum  number  possible,  so  the  amount  of  com¬ 
putational  work  required  is  bounded. 

As  just  described,  the  waves  are  used  to  update 
cell  values  directly.  In  the  actual  implementation,  the 
waves  are  used  to  update  the  fluxes  at  the  cell  interfaces 
by  calculating  the  flux  through  each  interface  due  to  the 
wave.  This  makes  these  boundary  conditions  easier  to 
use  in  conjunction  with  an  arbitrary  flux  differencing 
method  away  from  the  boundaries.  The  flux  at  each 
interface  is  modified  by  any  wave  that  crosses  the  inter¬ 
face.  For  example,  the  reflected  wave  in  Figure  3b 
crosses  four  interfaces  and  would  modify  the  flux  at  each 
of  these  interfaces. 

Finally,  we  must  calculate  the  flux  Hiti  at  the  solid 
wall  boundary  itself.  The  basic  flux  is  computed  by  solv¬ 
ing  a  Riemann  problem  at  the  wall  with  data  given  by 
Ui,j  and  Ditj.  The  vector  0,j  agrees  with  Ui%j  in  density, 
pressure  and  tangential  velocity  but  again  has  the  normal 
velocity  negated.  Ibis  flux  is  then  modified  by  any  wave 
that  reflects  off  the  wall,  once  by  the  outgoing  flux  of  the 
wave  AU  and  then  again  by  the  incoming  flux  of  the 
reflected  wave. 
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Figure  3  (a)  The  wave  completely  sweeps  through 
the  small  boundary  cell,  (b)  The  wave  reflects  off 
the  boundary  back  into  the  domain. 

4,  Mesh  Refinement 

The  adaptive  mesh  refinement  algorithm  (hen¬ 
ceforth  AMR)  is  based  on  the  use  of  uniform,  local  grid 
refinements  superimposed  on  an  underlying  coarse  grid. 
These  embedded  grid  refinements  can  be  recursively 
nested  to  maintain  a  fixed  level  of  accuracy  in  the  calcu¬ 
lation.  Unlike  other  embedded  grid  refinement  methods, 
(e.g  (181),  in  this  method  the  grid  cells  requiring 
refinement  in  each  level  are  grouped  together  into  rec¬ 
tangular  blocks  which  rc  uniformly  refined.  This  means 
that  some  coarse  grid  cells  may  be  unnecessarily  refined, 
but  has  the  advantage  that  all  grids  are  uniform  and  rec¬ 
tangular.  This  allows  us  to  maintain  vectorization 
without  using  gather/sca:ter  operations.  It  also  allows 
for  a  simple  user  interface,  since  a  finite  difference 
scheme  can  be  written  for  a  uniform  rectangular  grid 
without  concern  for  the  connectivity  of  each  cell.  The 
use  of  fine  grids  instead  of  unstructured  grid  points  also 
reduces  the  storage  overhead,  which  is  on  a  per  grid 
basis  for  our  method,  rather  than  the  overhead  per  cell 
found  in  unstructured  mesh  calculations.  The  additional 
complications  introduced  by  this  approach  occur  at  the 
interfaces  of  the  fine  and  coarse  grids  (see  below). 
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In  addition  to  refining  the  spatial  grid  for  time- 
accurate  computations  we  use  a  smaller  time  step  on  the 
fine  grids  as  well.  This  keeps  the  mesh  ratio  of  time  step 
to  space  step  the  same  on  all  grids,  and  so  the  same 
explicit  finite  difference  scheme  is  stable  on  all  grids. 
The  computational  work  is  thus  further  concentrated  on 
the  fine  grids,  where  it  should  be.  In  contrast,  some 
adaptive  methods  for  transient  flows  use  the  same  lime 
step  for  the  whole  mesh  [16,17],  This  can  be  less 
efficient,  since  the  resulting  Courant  number  may  be  far 
smaller  than  necessary  over  the  unrefined  portion  of  the 
grid. 

AMR  uses  an  automatic  error  estimation  pro¬ 
cedure,  based  on  Richardson  extrapolation,  to  determine 
the  regions  in  the  domain  where  the  resolution  in  the 
solution  is  insufficient.  These  coarse  grid  cells  are 
"flagged"  as  needing  refinement  In  addition,  the  irregu¬ 
lar  grid  cells  at  solid  bodies  should  be  flagged  as  needing 
refinement  if  the  geometry  of  the  boundary  is  under¬ 
resolved.  An  automatic  grid  generation  algorithm 
groups  these  flagged  cells  into  rectangular  grid  patches. 
We  have  developed  heuristic  procedures  that  are  quite 
successful  at  this  type  of  grid  generation  [4],  We  try  to 
balance  the  conflicting  goals  of  minimizing  the  number 
of  fine  grids  and  minimizing  the  area  that  is  unneces¬ 
sarily  refined. 

The  time  accurate  integration  algorithm  proceeds 
by  taking  one  step  on  the  coarsest  grid,  and  as  many 
steps  as  necessary  on  the  finer  level  grids  until  they 
catch  up  to  the  coarse  grid  time.  If  there  are  several  lev¬ 
els  of  fine  grids,  this  is  applied  recursively.  At  this  point 
the  grids  are  advanced  independently  of  each  other, 
except  that  fine  grids  require  boundary  values  from  adja¬ 
cent  fine  grids  or  interpolated  from  the  coarser  grids. 
For  a  five  point  stencil,  a  fine  grid  will  need  2  points  all 
around  the  outside  of  the  grid  in  order  to  advance  the 
solution  one  step.  If  there  is  an  adjacent  fine  grid,  it  can 
supply  the  missing  points.  Otherwise,  these  so-called 
dummy  points  are  obtained  using  bilinear  interpolation 
in  space  and  linear  interpolation  in  time  from  the  coarse 
grid. 

Since  we  will  be  computing  discontinuous  solu¬ 
tions  of  hyperbolic  conservation  laws,  the  adaptive  mesh 
refinement  algorithm  needs  to  be  conservative.  This  is 
complicated  by  the  use  of  different  time  steps  on  the  dif¬ 
ferent  grids.  Conservation  is  ensured  in  three  different 
pans  of  the  mesh  refinement  algorithm.  When  two  adja¬ 
cent  levels  of  grids  are  at  same  time,  the  fine  grid 
updates  the  coarse  grid,  performing  the  conservative 
averaging  procedure 
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where  r  is  the  mesh  refinement  ratio,  for  each  coarse  cell 


(i.j)  containing  a  fine  grid  cell  (*,/)  in  the  lower  left 
comer.  If  a  fine  grid  is  then  removed,  the  total  mass  in 
the  domain  is  conserved.  Secondly,  after  every  integra¬ 
tion  step  the  solution  is  post-processed  at  all  coarse  grid 
points  adjacent  to  a  fine  grid.  The  initial  coarse  flux 
(computed  ignoring  the  fine  grids)  is  subtracted,  and  the 
sum  of  the  fine  grid  fluxes  over  space  and  time  is  added 
in  its  place.  Thirdly,  we  use  conservative  interpolation 
procedures  to  initialize  the  solution  when  a  fine  grid  is 
created.  A  more  complete  description  oi  the  algorithm 
for  time-dependent  pdes  is  in  [6]. 

5.  Numerical  Results 

We  illustrate  the  method  on  two  time-dependent 
problems  involving  shock  waves.  In  the  first  example 
we  compute  flow  around  two  cylinders.  An  incident 
shock  travels  at  Mach  number  2.81.  One  cylinder  is 
slightly  ahead  of  the  other.  This  leads  to  an  interesting 
pattern  of  wave  reflections  between  the  two  cylinders.  In 
addition  there  is  a  reflected  bow  shock,  and  complicated 
wave  structures  behind  the  cylinders  after  the  shocks 
pass  by.  All  of  these  regions  use  the  adaptive  refinement 
as  the  solution  develops.  Figure  4  shows  the  incident 
shock  with  the  location  of  the  refined  grids  indicated. 
The  initial  coarse  grid  is  64  by  64.  Two  levels  of  grids 
are  used,  with  a  refinement  factor  of  4.  Figure  S  shows 
density  contours  of  the  solution  at  later  stages  of  the 
simulation. 
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Figure  4  Density  contours  and  grid  locations  at  in¬ 
itial  time  for  shock  impinging  on  two  cylinders. 

The  second  example  we  consider  is  a  Mach  2.2 
shock  travelling  in  a  channel  with  a  90  degree  bend. 
This  has  previously  been  studied  in  [1,21],  We  use  a 
very  coarse  initial  grid,  since  much  of  the  initial  tec- 
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Figure  5  Density  contours  at  later  times. 

(angular  grid  is  outside  the  computational  domain.  We 
use  two  additional  levels  of  refinement  by  a  factor  of  4  in 
each  case  in  order  to  obtain  good  resolution  of  the  shock 
and  induced  wave  pattern.  Figure  6  shows  density  con¬ 
tours  when  the  shock  has  passed  most  of  the  way 
through  the  channel.  In  this  figure,  the  location  of  the 
three  levels  of  refined  grids  is  indicated  on  the  contour 
plots. 

6.  Conclusions 

This  work  demonstrates  the  feasibility  of  an  adap¬ 
tive  Cartesian  grid  approach  for  fluid  problems  in  com¬ 
plicated  geometries.  Several  aspects  of  this  approach  are 
still  under  development  We  would  like  to  improve  the 
boundary  scheme  to  make  it  second  order  along  with  the 


Figure  6  Density  contours  for  shock  in  a  channel 

with  90  degree  bend. 

interior  scheme.  We  also  hope  to  eliminate  the  smooth¬ 
ing  step  now  used  to  insure  stability  in  the  very  small 
cells. 

We  also  plan  to  include  an  acceleration  procedure 
for  steady  state  calculations.  For  nested  Cartesian  grids 
multigrid  is  particularly  attractive,  since  the  data  struc¬ 
tures  and  grid  transfer  operations  in  the  adaptive  grid 
refinement  algorithm  make  up  almost  all  of  those  needed 
in  multigrid  [5].  Local  grid  refinement  and  multigrid 
have  already  been  combined  using  a  logically  rectangu¬ 
lar  body-fitted  grid  in  [7], 

An  important  consideration  is  whether  these  tech¬ 
niques  will  extend  to  Navier-Stokes  calculations.  For 
problems  with  boundary  layers,  it  is  usually  desirable  to 
use  body-fitted  grids  and  refine  heavily  in  the  direction 
normal  to  the  boundary.  Refining  Cartesian  grid  cells 
near  such  a  boundary  may  be  highly  inefficient.  In  such 
cases  a  component  grid  approach  may  be  useful,  in 
which  there  are  several  grids  with  distinct  coordinate 
systems.  For  example,  there  may  be  a  thin  body-fiued 
boundary  layer  region  in  addition  to  an  underlying  Carte¬ 
sian  grid.  These  multiple  components  will  overlap  in  an 
arbitrary  way,  creating  small  irregular  cells  as  in  the 
Cartesian  mesh  method  above.  Again,  stable  and  conser¬ 
vative  difference  equations  are  needed  to  compute  the 
flow  at  these  mesh  junctions.  The  techniques  used  here 
should  be  directly  applicable  to  this  situation.  Such  an 
approach  has  previously  been  considered  by  others,  e.g. 
[2,3).  However,  these  previous  efforts  did  not  treat  the 
interface  conditions  between  the  different  grids  in  an 
accurate  or  conservative  way. 
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